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COMPUTATION OF COMPLEX 
THREE-DIMENSIONAL TURBULENT FREE JETS 

Robert V. Wilson' and Ayodeji O. Demuren 2 
Department of Mechanical Engineering 
Old Dominion University, Norfolk, VA 23529 

ABSTRACT 

Three -dimensional, incompressible turbulent jets with rectangular and elliptical 
cross-sections are simulted with a finite-difference numerical method. The full Navier- 
Stokes equations are solved at low Reynolds numbers, whereas at high Reynolds numbers 
filtered forms of the equations are solved along with a sub-grid scale model to approximate 
the effects of the unresolved scales. A 2-N storage, third-order Runge-Kutta scheme is 
used for temporary discretization and a fourth-order compact scheme is used for spatial 
discretization. Although such methods are widely used in the simulation of compressible 
flows, the lack of an evolution equation for pressure or density presents particular difficulty 
in incompressible flows. The pressure-velocity coupling must be established indirectly. It 
is achieved, in this study, through a Poisson equation which is solved by a compact 
scheme of the same order of accuracy. The numerical formulation is validated and the 
dispersion and dissipation errors are documented by the solution of a wide range of 
benchmark problems. Three-dimensional computations are performed for different inlet 
conditions which model the naturally developing and forced jets. The experimentally 
observed phenomenon of axis-switching is captured in the numerical simulation, and it is 
confirmed through flow visualization that this is based on self-induction of the vorticity 
field. Statistical quantities such as mean velocity, mean pressure, two-point velocity spatial 
correlations and Reynolds stresses are presented. Detailed budgets of the mean momentun 
and Reynolds stresses are presented. Detailed budgets of the mean momentum and 
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Reynolds stress equations are presented to aid in the turbulence modeling of complex jets. 
Simulations of circular jets are used to quantify the effect of the non-uniform curvature of 
the non-circular jets. 
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Chapter 1 


INTRODUCTION 

1.1 Motivation 

Turbulent jets are present in many physical processes and technological applications. 
Turbulent jets can be found in combustors where the fuel and oxidizer are introduced as 
co-flowing jets, where the efficiency of such a process is largely determined by the mixing 
of the jets. Recently, jet aircraft noise has received much attention due to plans for a high- 
speed civil transport. A critical issue for the project’s success is reducing jet noise to 
acceptable levels near populated areas. The belief is that acoustic patterns can be altered 
by manipulating the large scale structures in turbulent jet flows through external forcing. 
Non-circular jets can also be used to enhance the mixing of hot jet gases with the sur- 
roundings in aerospace applications and thus avoid aircraft detection. In industrial applica- 
tions, efficient mixing is required to mix pollution issuing from smokestacks with the 
ambient surroundings to avoid its harmful effects. 

In the laboratory, turbulent jets usually originate from a high pressure stagnation 
chamber. Typically, the flow is then expanded through either a contoured nozzle or an ori- 
fice plate which caps the stagnation chamber. The jet is then allowed to mix with the ambi- 
ent surroundings and to develop in the streamwise direction. 

Experiments have shown that three-dimensional (3-D) jets can be used to enhance 
mixing and entrainment rates compared to nominally two-dimensional (2-D) jets. A 
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fundamental understanding of the dynamics of complex, turbulent jets is required for their 
prediction and control. The present study is concerned with understanding the spatial 
evolution of incompressible 3-D jets in the near to medium field. 


1.2 Survey of Previous Work - Experimental 

Early experimental studies of three-dimensional, turbulent jets issuing from nozzles 
and orifices (Sforza et al. 1966, Trentacoste and Sforza 1967, and Sfeir 1976) revealed a 
phenomenon known as axis switching whereby the orientation of the jet major and minor 
axes at the nozzle exit switch at a downstream location. 

Sforza et al. (1966) and Trentacoste and Sforza (1967) studied the mean flow of 3-D 
jets issuing from round, elliptical, and rectangular orifices of various aspect ratios. They 
characterized the streamwise development of the mean velocity using three distinct 
regions: (i) potential core, (ii) characteristic decay, and (iii) axisymmetric decay regions. 
In the potential core region, the mixing layer separating the jet core from the ambient 
surroundings at the orifice exit, has not spread to the jet centerline. As a result, the 
streamwise velocity near the jet centerline is constant in this region. In the second region, 
the velocity profiles in the plane containing the minor dimension of the orifice were found 
to be similar whereas those in the major plane are non-similar. Because the decay of 
centerline velocity was found to be dependent on orifice geometry, this region is referred 
to as the characteristic decay region. A third region is characterized by an axisymmetric 
decay of the centerline velocity which is proportional to the inverse of the streamwise 
coordinate. Velocity profiles in both major and minor planes were found to be similar and 
mostly independent of initial geometry. The results show that the length of the potential 
core region is roughly 5 diameters for rectangular and elliptical geometries of aspect ratio, 
AR = 10. The start of the axisymmetric decay region was roughly 50 diameters 
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downstream of the exit. Saddle-shaped streamwise velocity profiles (the maximum value 
occurs away from the jet centerline) were observed in the minor axis plane in the 
characteristic decay region. One axis switching was reported at 40 equivalent diameters 
downstream of the exit for the rectangular orifice of AR =10. 

Sfeir (1976) extended the earlier results by studying rectangular jets issuing from both 
nozzles and orifices for aspect ratios of AR = 10, 20, and 30. Axis switching and saddle- 
shaped velocity profiles were found for both orifice and nozzle jets. In general, jets from 
orifices showed axis switching locations closer to the exit when compared to jets from 
nozzles of equal aspect ratio. The saddle-shaped velocity profiles were more pronounced 
for jets from orifices. 

Krothapalli et al. (1981) performed experiments of rectangular jets from a moderate 
aspect ratio nozzle, AR = 16.7 at Reynolds number of 12,000. In the characteristic or two- 
dimensional region, they found self-similar profiles for the mean velocity, Reynolds shear 
and normal stresses and a linear growth of the jet width in the minor axis plane. The shape 
of the self-similar profiles was found to be insensitive to aspect ratio, for AR > 10. How- 
ever, the location where the self-similar profiles begin was found to be directly influenced 
by aspect ratio. Non-similar profiles were found in the major axis plane. 

Tsuchiya et al. (1985) studied the effect of exit shape on the mean velocity field of 
rectangular jets of aspect ratios, AR = 2 and 5. Smoothly contoured nozzles of various 
lengths and sharp-edged orifices were utilized as exit shapes. Only the jets issuing from 
orifice configurations produced saddle-shaped velocity profiles. All three configurations 
produced at least one axis switching event with the orifice jet the closest to the jet exit. In 
a later study, Tsuchiya et al. (1989) reported the saddle-shaped velocity profiles in nozzle 

jets as well. 
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Ho and Gutmark (1987) studied jets issuing from elliptic nozzles with AR =2. Entrain- 
ment rates for the elliptic jet were found to be several times larger than equivalent area 
plane or axisymmetric jets. The increased entrainment rate was explained in terms of vor- 
tex induction due to the non-uniform azimuthal curvature of the shear layer. Fluid was 
found to be preferentially entrained at the minor axis plane as this portion of the vortex 
moved outward, thus resulting in an axis switching. Three such events were observed 
before the jet approached a circular shape. 

Hussain and Husain (1989) extended the study of elliptic jets to include the effect of 
initial condition of the boundary layer at the jet exit plane, such as azimuthal variation of 
momentum thickness, turbulence level and effect of forcing. Axis switching was reported 
for up to 100 equivalent diameters. Subsequent studies explained the vortex ring pairing 
process (Husain and Hussain 1991) and the preferred mode coherent structure (Husain and 
Hussain 1993). 


Zaman (1996) used azimuthal and streamwise vorticity dynamics to explain the pres- 
ence (or absence) of axis-switching in low aspect ratio, AR =3, rectangular jets. The study 
also investigated the effect of adding vorticity generating tabs at the nozzle exit. It was 
shown that contracting nozzles upstream of the jet exit plane could produce two pair of 
counter-rotating streamwise vortices which eject fluid form the jet core to the ambient. 
This sense of rotation did not promote axis-switching within the measurement domain. 
Tabs placed on the short sides of the rectangular jet produced two pair of streamwise vor- 
tices which pump fluid from the ambient to the jet core resulting in rapid axis switching. 
Tabs placed on the long sides of the jet produced streamwise vortices of the same sense as 
that from contracting nozzles, resulting in no axis-switching. 
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Periodic forcing of the jet was shown to spatially organize the azimuthal vortex struc- 
tures resulting in the first axis-switching location being closer to the jet exit when com- 
pared to the naturally developing jet. This was used to explain axis switching in supersonic 
screeching jets. The authors note that the effects of azimuthal and streamwise vorticity are 
not mutually exclusive and that the effect of streamwise vorticity pairs can either stop or 
augment axis-switching. 

1.3 Survey of Previous Work - Analytical 

Analytical techniques have been used to study non-circular jets in the literature. These 
methods fall roughly into two categories: (i) linear stability analysis, and (ii) vortex meth- 
ods. In the former method, the governing equations are linearized about a base flow and 
the perturbation quantity is assumed to be composed of normal modes. The linear stability 
analysis of non-circular layers is more complex than nominally 2-D, plane shear layers 
due to the inherent three-dimensionality of the base flow. 

The stability of elliptic jets was initially studied by Crighton (1973) using a vortex 
sheet model. Morris (1988) extended the analysis to finite thickness shear layers. 
Koshigoe and coworkers (Koshigoe and Tubis 1986, 1987, Koshigoe et al. 1987) studied 
the instability of circular and elliptical jets using a generalized shooting method. The stud- 
ies showed that for elliptic jets, instabilities associated with the lower curvature portion of 
the jet boundary layer are dominate inferring that large scale coherent structures would 
form first in the minor axis plane. 

Tam and Thies (1993) investigated instability waves of rectangular jets using a vortex 
sheet model which approximates the region very close to the jet exit where the boundary 
layer thickness is very small. The analysis identified four linearly independent families of 
instability modes based on mode shape (symmetry considerations). The authors found that 
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within each family, the first and third modes are associated with jet comer instability. The 
second mode is associated with instability of the mixing layers of the center of the flat 
sides of the jet. The second mode was found to have the largest spatial growth rate and is 
thus expected to be the dominate instability in rectangular jets. 

The dynamics of isolated inviscid elliptic vortex rings was studied by Viets and Sforza 
(1972) and Dhanak and DeBemandis (1981) in an effort to model elliptic jets from noz- 
zles. An analysis using the Biot-Savart law predicts induced velocity perpendicular to the 
plane containing the vortex core. The magnitude is proportional to the local curvature of 
the vortex tube, C, and the log of the inverse of the cross-sectional area, A\ it - c£log(/C‘) , 
where t> is the local bi-normal vector. Therefore, vortex tubes with high local curvature 
and small cross-sectional areas will experience larger induced velocities. The analysis pre- 
dicts that an initially planar elliptic ring will become distorted and switch its major and 
minor axis. Axis switching in jet issuing from nozzles and orifices is considerably more 
complex due to viscous and turbulent diffusion, shear, entrainment, and flow instabilities. 

1.4 Survey of Previous Work - Computational 

Previous computational studies of three-dimensional free jets are reviewed in this sec- 
tion. In comparison to the available experimental studies, the number of computational 
studies in the literature is sparse. 

Early attempts at a numerical solution of 3-D jet flow utilized the Reynolds Averaged 
Navier-Stokes (RANS) equations whereby the instantaneous Navier-Stokes equations are 
first time-averaged and a turbulence model is used to close the system of equations (the 
well-known closure problem). The RANS equations are then solved for the time-averaged 
velocity and pressure. This procedure was followed by McGuirk and Rodi (1979) in their 
study of 3-D free jets of aspect ratio, AR = 1,5, 10, and 20. A k-e turbulence model was 
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used to close the system of equations. The flow was assumed to be parabolic in the stream- 
wise direction, allowing the use of a spatial marching procedure in the downstream direc- 
tion. While this procedure requires less CPU time and memory than an elliptic solution, a 
boundary layer analysis shows that the flow is only parabolic in the far field. The computa- 
tions were unable to reproduce the experimentally observed axis switching and saddle- 
shaped velocity profiles. One axis switching event was predicted only after an ad-hoc 
specification of the lateral velocity components at the inflow of the computational domain 
was made. The saddle-shaped velocity profiles were never predicted. The authors attrib- 
uted this deficiency to the Jfc-e turbulence model which is incapable of capturing the 
effects of turbulence driven secondary motion. 

Quinn and Militzer (1988) utilized a 3-D elliptic solution procedure to solve the 
RANS equations for the turbulent square jet from a sharp-edged orifice. Unlike the 
McGuirk and Rodi study, the velocity components at the inflow were specified from the 
author’s experimental results, which were also presented in the paper. The computations 
were successful in predicting the decay of the centerline velocity in the medium to far 
field. Results in the near field, x/D e < 5, were only qualitatively predicted, which the 
authors attributed to a relatively coarse grid. The experimental results showed off-center- 
line peaks of mean streamwise velocity in the near field, a faster spread rate when com- 
pared to a circular jet of equivalent area, and positive mean static pressure in the very near 
field at the jet centerline. Saddle-shaped profiles were also reported for the normal Rey- 
nolds stress components. Detailed velocity profiles from the numerical solution were not 
provided so that a comparison with the above mentioned trends is not possible. 

Only in the last decade an unsteady numerical solution of the Navier-Stokes equations 
for 3-D jet flows has been possible. In the direct numerical simulation (DNS) approach, all 
scales of motion are resolved by the computational grid and no modeling is required. In 
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the large eddy simulation (LES) approach, the energy containing large scales of motion 
are resolved, while the unresolved small scales are modeled. 

Grinstein and DeVore (1992) performed LES of spatially-developing square jets at 
moderate Reynolds numbers to study large scale coherent structures. The Euler equations 
of motion were solved (i.e. molecular viscosity is neglected) and an explicit filtering of the 
velocity was used as a minimal subgrid model for the unresolved scales of motion. The 
authors explain the jet dynamics in terms of the deformation, merging, and breakdown of 
initially planar square vortex rings. The deformation of the initially planar rings was 
attributed to the induced velocity due to azimuthal curvature present in the jet comers. The 
relative induced velocity results in the comers of the square ring moving ahead and 
towards the jet centerline, while the flat sides remain behind and move away from the cen- 
terline. The deformation results in a switching of the orientation of the square jet by 45° at 
an axial location downstream of the jet exit (x/D e ~ 0.8- 1.0). Flow visualization of the 
results revealed pairs of counter-rotating streamwise hairpin vortices in the high strain cor- 
ner region between two adjacent vortex rings. Pairing of the vortex rings was accompanied 
by amalgamation of neighboring hairpin vortices which doubled their streamwise extent 
and led to the eventual breakdown of the rings. Subsequently, the flow was characterized 
by less organized, small scale vortices, indicative of fully turbulent flow. 

Later studies (Grinstein 1993) focused on the LES of the very near field ( x/D e < 5) of a 
2: 1 aspect ratio rectangular jet. Thus, the experimentally observed axis switching at x/D e = 
7 is not captured in the author’s computation. Other studies focused on the vorticity 
dynamics of isolated, rectangular vortex rings (Grinstein 1995) and the effects of com- 
pressibility and initial condition (Grinstein 1996). 
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Miller et al. (1996) performed simulations of non-circular jets at low Reynolds number 
(Rep = 800) for a streamwise extent of x/D e = 9. The jets are forced at the inflow plane 
with a single sinusoidal mode. Iso-surfaces of instantaneous vorticity magnitude show the 
flow to be laminar - i.e. composed of smooth, symmetrical structures. Axis-switching is 
predicted at x/D e = 3.1 for the 2:1 AR elliptical jet and x/D e = 6.3 for the 2: 1 AR rectangu- 
lar jet. Plots of centerline velocity reveal that the end of the potential core is not reached 
within the computational domain (0-9 D e ). 

1.5 Summary of Literature Survey 

General conclusions from the experimental and computational studies on non-circular 
jets in the literature are made in this section. Even though there are some conflicting opin- 
ions, trends are consistently observed in many experiments. 

Early studies showed that the decay of centerline velocity in rectangular jets was char- 
acterized by three distinct regions: (i) potential core region, ending at roughly four to five 
diameters, where the centerline velocity is constant, (ii) characteristic decay region, end- 
ing at roughly 20 - 60 diameters (dependent on aspect ratio), where decay is dependent on 
initial geometry and profiles in the minor axis plane only are similar, and (Hi) the axisym- 
metric region where the decay is proportional to the inverse of the streamwise coordinate 
and is mostly independent of initial condition. There is evidence that for jet with AR > 10, 
the axis switching location scales linearly with nozzle aspect ratio. The strong skewing of 
streamlines near the jet exit in orifice jets results in axis switching closer to the jet exit in 
comparison with equivalent nozzle jets. Pronounced saddle-shaped profiles for streamwise 
mean and fluctuating velocity were also observed with orifice jets. Less pronounced sad- 
dle-shaped profiles were observed for nozzle jets as well. Studies reveal that entrainment 
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and mixing in non-circular jets can be several times larger in comparison with equivalent 
area, circular jets. This increase occurs preferentially in the minor axis plane. 

Linear stability analysis of rectangular jets indicate that instabilities associated with 
the center of the flat sides of the jet have larger spatial growth rates compared with comer 
modes and are thus expected to be dominant. 

Numerical studies at lower Reynolds numbers (Re D ~ 800 ) where the flow is laminar/ 
transitional, predict axis switching of non-circular jets forced with a single sinusoidal 
mode. The effect of the spectral content of the forcing function has not been addressed in 
such studies. Numerical simulation at higher Reynolds numbers are required for compari- 
son with most jet experiments. Spatially-developing LES performed at higher Reynolds 
numbers are limited to the potential core region and single sinusoidal mode forcing. Axis 
switching for the rectangular or elliptic jet has not been simulated numerically at higher 
Reynolds number. 


1 .6 Objectives of the Current Study 

The specific objectives of the current study are given in this section which attempt to 
address issues not covered in the literature. The objectives of the current study are enumer- 
ated below. 

(i) Develop a higher-order accurate numerical formulation for the simulation of spatially- 
developing, unsteady, incompressible flows with improved resolution of high fre- 
quency modes. 

(ii) Show the effect of initial condition on jet dynamics by altering the spectral content of 
the forcing function. 
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(Hi) Simulate the potential core and characteristic decay regions of non-circular jets at 
higher Reynolds numbers (Re D - 10 5 ). 

(iv) Demonstrate axis switching at higher Reynolds numbers and explore the axis switch- 
ing mechanism. 

(v) Compute budget terms of the mean momentum and Reynolds stress equations which 
can be used in the turbulence modeling of complex jets. 

A description of the organization of this study is now given. The governing equations 
for incompressible flow are presented in Chap. 2 along with the necessary subgrid stress 
model required for large eddy simulation. Boundary and initial conditions for the spa- 
tially-developing jet are also presented. The temporal and spatial discretization of the gov- 
erning equations are presented in Chap. 3, while the details of the solution of the Poisson 
equation for pressure are presented in Chap. 4. The numerical formulation is validated in 
Chap. 5 through the solution of a variety of benchmark problems. Results of the direct 
numerical simulation of rectangular jets are presented in Chap. 6. In Chap. 7, results from 
the large eddy simulation of rectangular, elliptic, and circular jets at higher Reynolds num- 
ber are presented. Finally, a summary and conclusions from the study are provided in 
Chap. 8. 
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Chapter 2 


MATHEMATICAL FORMULATION 

The equations governing the conservation of mass and momentum of an isothermal, 
incompressible, time-dependent fluid are developed in this section. The initial and bound- 
ary conditions for the spatially-developing jet are also given. 

2.1 Governing Equations 

Figure 2. 1 shows the computational domain for the jet simulations along with the 
coordinate system and domain dimensions. The inflow boundary of the computational 
domain is located at a finite distance downstream of a hypothetical nozzle which produces 
a thin boundary layer separating the jet core from a stagnant freestream. It is assumed that 
the streamwise velocity makes a smooth transition from the jet velocity at the core to the 
ambient velocity and is thus modeled with the hyperbolic tangent function. The fluid then 
leaves the computational domain through the outflow boundary located a distance, L x , 
downstream of the inflow boundary. 

Freestream boundaries are located in the y and z coordinate directions where fluid is 
entrained into the jet. Unless otherwise specified, the governing equations and reported 
quantities will be normalized using the equivalent jet diameter at the inflow plane, D e , as 
the length scale and the jet core velocity at the inflow plane, U a , as the velocity scale. The 
equivalent diameter for the non-circular jet is defined as the diameter of a circle having the 
same area at the inflow plane. Time will be normalized using the time scale, D/IJ 0 . 
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2.1.1 Continuity and Momentum Equations in Physical Space 


The Navier-Stokes equations written in non-conservative form for an incompressible 
fluid are given in this section. The continuity equation in a Cartesian coordinate system 
written in indicial notation is given by: 


where Xj, uj,j = 1 , 2 , 3, represent the spatial coordinates and instantaneous velocities in the 
physical Cartesian coordinate system. 

The non-dimensionalized momentum equation in a Cartesian coordinate system is 
given by: 


d«, _ _dp 1 d“i 

dt +U >dXj dx i ; + Re D dxjdxj 


(2.2) 


where Re D = ( U a D e )/\ is the flow Reynolds number, p is the non-dimensional pressure, 
and v the kinematic viscosity of the fluid. The momentum equation has three scalar com- 
ponents (i = 1,2,3). 


The large eddy simulation (LES) approach is explored in this study, in which the large 
scales of turbulent motion are resolved while the smallest scales are not computed directly 
and are modeled in terms of the resolved scales. The filtering operation is defined by: 


/(*) = J/($)G(i -$)</$ (2.3) 

where the integration is extended over the computational domain, ft, and the general vari- 
able /is filtered to yield the spatially averaged value, f. The variable G, denotes a spatial 
filter which must satisfy the normalization constraint: 
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( 2 . 4 ) 


Jg<* -$)</$ - i 

n 

It is convenient to define a grid filter that commutes with spatial and temporal differen- 
tiation, such that: 


u = u 

dt dt 


( 2 . 5 ) 


V = U 

dx dx 


( 2 . 6 ) 


Applying the grid filtering operation of Eq. (2.3) to the governing equations of motion, 
Eqs. (2.1) and (2.2), the filtered equations of motion are obtained as 


duj 

d Xj 


= 0 


( 2 . 7 ) 


+ = dQ ! 1 dfi, dq tj 

dt J dxj dx, Re D dx ) dx J dx J 


( 2 . 8 ) 


where $ = p + (\/3)x tk is the pseudo pressure, = u i u j -u i is the unresolved subgrid 
scale stress due to the non-linearity of the convection terms, and q Lj = x i . - (l/3)8 ;> t M is the 
anisotropic part of the subgrid scale stress. The subgrid scale (SGS) models considered in 
this study are defined in the next section. 


2.1.2 Subgrid Scale (SGS) Models 


The purely dissipative model of Smagorinsky (1963) is used in the current study. The 
purpose of the SGS model is to account for the unresolved small scales. The Smagorinsky 
model has been applied to the LES of many turbulent flows such as homogenous isotropic 
flow, channel flow, and mixing layer flow. The Smagorinsky model has been one of the 
most popular SGS models for LES, partly because it correctly models the global transfer 
of energy from large to small scales. It provides an energy sink such that the large scale 
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energy is drained from the flow. However, it does not correctly model local effects such as 
solid boundaries and near-wall regions, localized transfer of energy from small scales to 
large scales, and laminar and transitional flows. 

The anisotropic part of the subgrid stress is modeled in terms of the resolved scales 
using an eddy viscosity model due to Smagorinsky (1963): 

«,j = -2v,S/> (2-9) 

where v, is the turbulent viscosity predicted by the SGS and the strain rate tensor of the 
resolved scales, given by 


5 * = 2 


^du [ du j 


The turbulent viscosity is given by 


( 2 . 10 ) 


v, = CA%| (2.11) 

where C = 0.01 , A 2 is the volume of the computational cell, and is the mag- 

nitude of the resolved strain rate tensor. 

2.1.3 Continuity and Momentum Equations in a Mapped Coordinate System 

The governing equations given in Sec. 2.1.1 are mapped to an alternate coordinate sys- 
tem through the use of the chain rule which introduces “metric” terms. This approach has 
the potential for a more efficient use of grid points in resolving the thin boundary layer at 
the domain inflow. The velocity components are defined using the Cartesian coordinate 
system while the spatial gradients are defined in terms of the computational coordinate 
system with uniform grid spacing. The gradient terms can then be discretized using high- 
order compact finite difference schemes. 
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Consider the following general mapping of Cartesian coordinates (x„jt 2 ,jr 3 ) to the 
alternate coordinates (e,, e 2 , e 3 ) : 

E m = E m(* I. *2. *3> (2-12) 

First derivatives in the Cartesian coordinate system are expressed in terms of deriva- 
tives in the alternate coordinate system using the chain rule: 


JL = 

dxj a* ; ae m 


( 2 . 13 ) 


Upon expansion, Eq. (2.13) represents a sum of three products of metrics and deriva- 


tives. The Laplacian operator in the alternate coordinate system is expressed as 


a _ ae m a fa e „ d ) 

dxjhxj ax y ae m [ax y ae„J (214) 

Equation (2.14) involves 18 terms upon expansion. It is understood that the term 
a/ae m operates on the term in parenthesis. Equations (2.13) and (2.14) are used to express 
spatial gradients appearing in the continuity and momentum equations in terms of gradi- 
ents in the alternate coordinate system. The continuity equation in the alternate coordinate 
system is given by 


ae_a« : 

« J_ 

dx: dt m 

J m 


= 0 


(2.15) 


The momentum equation in the alternate coordinate system is given by 


du, t u de^du, _ _<Kdp_ + J 

dt Ui dxjdz m 3x ( 3e m + Re D dxjdE m \d Xj dE nj 
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The non-linear terms present in the momentum equation preclude the exact solution of 
the governing equations which must be solved numerically. The numerical approximations 
to the governing equations are given in Chap. 3. 

The continuity and momentum equations represent four scalar equations for the four 
unknown variables (three velocity components in the Cartesian coordinate system and 
pressure). In the next chapter it is shown that the continuity equation is enforced through 
the solution of a Poisson equation for pressure. The Poisson equation is obtained by taking 
the numerical divergence of the discretized momentum equation. 

2.2 Boundary Conditions 

In order to define a well posed problem, the boundary and initial conditions for the jet 
simulations are defined. The ellipticity in the spatial terms of the governing equations 
requires that boundary conditions be defined at all boundaries. A diagram of the boundary 
conditions is provided in Fig 2. 1 . 

2.2.1 Stream wise Inlet Boundary Conditions 

In the laboratory, jet flows are commonly generated by the use of a fan which forces 
fluid along an enclosed nozzle. The jet leaves the exit plane of the nozzle where it interacts 
with the ambient fluid. Prior to exit, the jet can be considered as a relatively uniform 
freestream and a curved boundary layer at the walls of the nozzle. A short distance down- 
stream of the nozzle exit, the boundary layer is smoothed so that the mean streamwise 
velocity can be modeled using the hyperbolic tangent (tanh) function. The inflow bound- 
ary of the computational domain is placed at a short distance downstream of the nozzle 
exit which is not actually included in the jet simulations. 
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(a) Mean Inlet Boundary Conditions 


The mean or time-averaged streamwise velocity component, U, at the inflow boundary 
is given by 


y<0.**>= «, + <^i- , U nh( | L) (2., 7) 

where U c = (U H + U L )/ 2 is the convective velocity and U H and U L represent the velocities 
of the jet core and ambient fluid, respectively. The quantity, r, represents the minimum 
directed distance from the point, ( 0 , y, z) to the line of constant convective velocity of the 
boundary layer (see Fig. 2.2). Thus, the shape of the line of constant convective velocity 
determines the jet geometry at the inflow plane. A super-elliptic equation was used to 
specify all jet geometries in this study: 



where a = AR*b and b represents the semi-major and minor dimensions of the constant 
convective velocity contour, and n is the exponent in super-elliptic coordinate system. For 
example, n = 2 defines an elliptic contour, while n >> 1 defines a rectangular contour with 
slightly rounded comers. The momentum thickness of the boundary layer at the inflow 
plane, 0„ is used to normalize the directed distance, r. If the point ( 0 , y, z) is “outside” the 
boundary layer contour as in Fig. 2.2, r is defined to be negative, while r is defined to be 
positive if the point lies on the inside of the boundary layer contour. Equation (2.17) pro- 
duces a constant thickness boundary layer if the momentum thickness, 0 O , is constant at 
all azimuthal positions along the boundary layer. Non constant thickness boundary layers 
are generated by specifying the desired variation of 0„ along the contour of the boundary 
layer. 
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The mean streamwise velocity at the inflow plane is specified by defining the jet core 
and ambient velocities (U H and U L ) and the momentum thickness of the shear layer, 9 0 in 
(2.17). For this study U H = 1, U L = 0.01, and D/Q 0 = 30, which models the experimental 
jet profile a short distance ( x/D e < 1) downstream from a contoured nozzle. For elliptic jets 
the value of the exponent of the super elliptic coordinate system is n = 2, while for the 
rectangular jet, n = 10 is used which slightly rounds the comers (see Fig. 7.9a). This also 
models experimental jets where viscous diffusion is known to smooth the shear layer exit- 
ing from rectangular nozzles at the sharp comers. The mean transverse velocity compo- 
nents at the inflow plane are set to zero. 

(b) Forced Inlet Boundary Conditions 

In an effort to model jet experiments, a time dependent forcing function of low inten- 
sity is added to the mean velocity components at the inflow boundary to promote unsteady 
motion. At higher Reynolds numbers and computational lengths, it is speculated that small 
round-off errors would grow to produce unsteady motion of the unstable shear layers, thus 
obviating the need for forcing functions. 

Two classes of perturbations are used in the current study; (i) sinusoidal perturbations 
of a specified frequency and (ii) perturbations having an experimentally measured velocity 
spectrum and transverse root mean square ( rms ) value, i.e., 

(i) Perturbations from linearized viscous stability theory 

The first class of perturbations is derived from the solution of the Orr-Sommerfeld 
equation (OSE) which governs the instability of the reference hyperbolic tangent profile to 
spatially developing disturbances. The details of the solution of the OSE are presented in 
Wilson and Demuren (1996) with the general form of the perturbation velocities being: 
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u(y,t) = Real{<() , (>-)e'< , >'}, v(y, t) = Real{-ia<t>(y)e““'} (2.19) 

where a , $(y) , to represent the perturbation wavenumber, eigenfunction, and angular fre- 
quency, respectively. The variables, y and t represent the transverse coordinate and time, 
respectively, while i = J -\ . The stability equations are solved for the plane mixing layer 
and the resulting eigenfunctions are adapted to non-circular jets by replacing the trans- 
verse coordinate, y, in Eq. (2.19) with the minimum directed distance, r. 

(ii) Random perturbations from experimental data 

Perturbations which have a broad spectrum are generated based on experimental data 
for a plane mixing layer. The velocity power spectra and root-mean-square ( rms ) perturba- 
tion levels were taken from the experiment of Spencer and Jones (1971). Because phase 
information is not included in the power spectra, a random phase relationship for the 
modes comprising the spectra was assumed. The velocity perturbations are found by per- 
forming a Fourier transform of the complex Fourier coefficients, dj(y,f), defined by: 

FAy, fjufT I 

o t (y,f) = J(cosy(y) + /sinY(y)) (2.20) 

where F t (y, f) , u} , T, and y(>) represent the normalized spectrum function, the rms level 
of the i th velocity component, time interval of the simulation, and the randomly generated 
phase angle, respectively. The velocity perturbations for the mixing layer represented by 
Eq. (2.20) are adapted to non-circular jets by replacing the transverse coordinate, y, with 
the minimum directed distance, r. The complete details of the derivation of time-depen- 
dent inlet boundary conditions based on a experimentally measured spectra are given in 
Wilson (1993). 
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The random inlet boundary conditions are the spatial analog to the random perturba- 
tions generated for initial conditions in temporal simulations. Figure 7.1a shows represen 
tative time traces and power spectra for the broad mode forcing function. 

It is widely accepted that the proper boundary condition for the Poisson equation for 
pressure (derived in the next chapter) is the Neumann boundary condition (Gresho and 
Sani 1987). This condition is derived by applying the normal component of the momen- 
tum equation at the boundary. Applying the i = 1 component of Eq. 2.2 and solving for the 
pressure gradient term results in Neumann condition for pressure at the inflow: 


[dpi _ 

du, 

3 u, ] 3 2 «i 

La*J. 

~di 

■ U] Jxj + Re De dxjdXj 


where the subscript “o” is used to denote the inflow boundary plane. The first term on the 
RHS of Eq. 2.21 is known from the inflow velocity boundary conditions, and the second 
and third terms are also known before the solution of the Poisson equation as will be 
shown in the next chapter. 

2.2.2 Streamwise Outflow Boundary Condition 

A characteristic analysis of the governing elliptic differential equations reveals no real 
characteristic curves along which disturbances travel. A disturbance is instead propagated 
in all directions at once. As a result, the solution of elliptic partial differential equations 
requires the specification of boundary conditions along the entire boundary. Boundary 
conditions are well defined at the inflow plane and can be reasonably approximated at the 
freestream boundary which is placed a large distance from the jet dynamics at the center- 
line. However, the conditions at the outflow boundary are not known a prior and must be 
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specified such that those errors do not adversely affect the accuracy in the interior of the 
domain. 

Sources of ell ipticity in the governing equations that link errors at the boundary to the 
interior can be traced to pressure and viscous terms. In addition, it is possible to convect 
errors from the outflow to the interior if flow reversal occurs. The source term for the Pois- 
son equation for pressure also links outflow errors to the interior. 

A technique which breaks the link between errors at the outflow boundary and the 
interior solution is the buffer domain technique proposed by Streett and Macaraeg (1990). 
This technique is adapted to the current jet simulations whereby the ellipticity of the gov- 
erning equations is eliminated in a small region (called the buffer layer) which is appended 
to the solution domain at the outflow. Practically, this is achieved by a factor which turns 
off those elliptic terms in the buffer layer. The resulting momentum equation becomes: 


dt 


+ [/««/+ (l-/«)£/ 7 )g^ 



fbl 

(<l N 

d U, 

1 


Re D \ 


+ Re Df 

U*2 dx ] , 


( 2 . 22 ) 


where fa is the buffer layer factor which gradually changes from unity in the solution 
domain to zero in the buffer layer through the following function: 


fbi(x l) = 1 - tanh[c M (jr, -x ]/2 )]} (2.23) 

where c bl = ln[/ cr /(l -/„)]/( l(x cr -x w2 )) is a constant which controls the rate of transi- 
tion between the solution and buffer domains while xy 2 controls the transition location. 
For the simulations presented in Chaps. 6 and 7, the following buffer layer parameters are 
used,/ cr = 0.99999, x cr = 0.99 L x , and x^ = 0.9L^. This results in a computational domain 
length of ten diameters and a buffer domain length of two diameters. The convection 
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velocity in the buffer domain. Up is computed from the time-averaged velocity as the com- 
putation progresses. In addition, the source term for the Poisson equation for pressure is 
gradually forced to zero in the buffer domain through the same function. Zero gradient 
boundary conditions for all velocity components and the pressure are then applied at the 
outflow. 

2.2.3 Freestream Boundary Conditions 

The freestream boundary of the jet simulations is located at a large distance from the 
jet centerline such that the developing structures do not cross this far-field boundary. This 
permits the freestream boundary to be modeled as a zero-traction boundary where pres- 
sure is set to zero and the gradients of the instantaneous velocity field are set to zero in the 
transverse and lateral directions: 

p = 0 ( 2 . 24 ) 



where n denotes the direction normal to the freestream boundary. 

2.3 Initial Conditions 

A spatial simulation of turbulent jet flow is performed in this study where a fixed 
region of the flow is computed and disturbances grow in the streamwise direction. This 
can be contrasted with a temporal simulation where a small region of the flow is followed 
in time and the domain moves in the streamwise direction. 

As a result of the spatial reference frame, initial conditions are of minor importance 
because they are quickly convected out of the domain and the dynamics of the jet flow are 
determined by the forcing functions applied at the inflow plane. Simulations are started on 
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coarse grids with the velocities specified at the inflow plane used to initialize the velocity 
field in the interior at t - 0. After several flow through times (the time required for a fluid 
particle to convect from the inflow to the outflow plane at the convective velocity, Uq) the 
initial conditions are “washed” from the domain. Simulations on finer grids are started 
from results on coarser grids by prolongating the results using a standard, second-order 
accurate interpolation formula. 
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Fig 2.1 Computational domain and coordinate system for the spatially-develop- 
ing jet (iso-surface of vorticity magnitude from LES of rectangular jet 
shown). 



Fig 2.2 Diagram showing the minimum directed distance, r, from the point 

(0, y, z ) to the contour of constant convective velocity, U c , at the inflow. 
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Chapter 3 


NUMERICAL FORMULATION 

The numerical approximations of the governing equations are described in this chap- 
ter. The temporal approximation to the governing equations is given in this section, while 
the spatial approximations are given in the next section. 

3.1 Temporal Discretization 

The time advancement scheme used to march the momentum equations in time should 
possess several qualities; low dispersion and dissipation errors over a wide range of step 
sizes, low-storage requirements, and a relatively large stability envelope. The family of 
low-storage, Runge-Kutta schemes proposed by Williamson (1980) possesses these desir- 
able qualities. The scheme is low-storage in the sense that only two storage locations (one 
for the time derivative and one for the variable itself) are required for time advancement. 
In comparison, a third-order fully implicit scheme requires four storage locations. For 
simplicity, the numerical approximations for the governing equations are given in the Car- 
tesian coordinate system with uniform grid spacing. Extension of the formulation to curvi- 
linear grids is accomplished by using the chain rule to replace the derivatives in physical 
space with derivatives in the uniform computational space and is straightforward. 

The additional metric terms are discretized using the same higher-order compact 
schemes. The momentum equation is advanced from time level, N, to JV+1, using Q 
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substages. The temporal derivative in the momentum equation is discretized using a third- 
or fourth-order explicit Runge-Kutta scheme: 


3 ,/f+l M ^ M 

~ ^ _ jjM op (2 I) 

dt ~ b M At 1 3 *- 

where b M is a constant of the Runge-Kutta scheme, uf represents the ith velocity compo- 
nent at the Mth substage (M = 0 is the Nth time level, M=Q is the N + 1 time level). The 
term, , denotes the sum of convection and diffusion terms: 


H 


M 


M f -.2 M \ 

du, 1 d_Uj 
Uj d*j + Re D r [ dx j dx j / 


( 3 . 2 ) 


The terms on the right hand side (RHS) of Eq. (3.1) are assumed known from the previous 
sub-stage or from the initial conditions at t = 0. The calculation of the pressure is accom- 
plished by solving a Poisson equation at each sub-stage such that the continuity is 
enforced. Since the pressure, p M , is calculated before the advancement of Eq. (3.1), «, , 

can be calculated explicitly using Eq. (3.1). 


The low-storage requirement is accomplished by continuously overwriting the storage 
location for the time derivatives and unknown variables at each sub-stage: 


fif <r- a M H*?- ' ( 3 . 3 ) 

u “+'^u“ + b M AlH, M ( 3 . 4 ) 

where fi“ = H* - [dp M /dx,] and the notation is used to indicate that the storage 

locations, H ?~' , are overwritten by, //**, «f +1 , respectively. Tables 3.1 and 3.2 show 
the constants, a M and b M (to 8 significant figures) for the low-storage, third- and fourth- 
order schemes, respectively. 
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Table 3.1 Coefficients of third-order Runge-Kutta schemes 
from Lowery and Reynolds (1986) 


M ‘ 

a" 

b M 

1 

0 

0.500 

2 

-0.68301270 

0.91068360 

3 

-1.33333333 

0.36602540 


Table 3.2 Coefficients of fourth-order Runge-Kutta schemes 
from Carpenter and Kennedy (1994) (defined there as solution 3) 


M 


b M 

1 

0 

0.14965902 

2 

-0.41789047 

0.37921031 

3 

-1.19215169 

0.82295502 

4 

-1.69778469 

0.69945045 

5 

-1.51418344 

0.15305724 


The stability characteristics of the Runge-Kutta schemes can be analyzed by consider- 
ing the model equation: 


gj *=£(♦,») (3.5) 

where $ is the generic unknown to be advanced in time and k is the time derivative which 
contains the spatial terms of the governing equation. Equation (3.5) is transformed from 
physical space to wavenumber space by decomposing <t> into Fourier modes: 

♦ = kt)e ikx (3.6) 

where <j>(r) is the Fourier coefficient of <j> , i = J-i , and k is the wavenumber. Substituting 
Eq. (3.6) into Eq. (3.5) yields: 
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(3.7) 


where X (a complex number) is the Fourier symbol of the spatial operator H . The Runge- 
Kutta scheme is used to expand the term on the LHS of Eq. (3.7) which gives an amplifica- 
tion factor, G = <j> A, + 1 /<j> N for the third-order scheme: 

G = 1 + (A. A/) + J (A»Af) 2 + 2(A.Ar) 3 (3.8) 

2 o 

It can be shown that all three-stage, third-order Runge-Kutta schemes have the same 
amplification factor given in Eq. (3.8). Solving Eq. (3.7) analytically in time gives the 
exact amplification factor, G e : 


G e = (3.9) 

Comparing Eqs. (3.8) and (3.9), the three-stage, third-order Runge-Kutta scheme is a 
polynomial approximation to the exact solution to third-order. Similarly, the five-stage, 
fourth-order Runge-Kutta scheme has amplification factor: 

G = 1 +(A.A/) + ^(AiAr) 2 +g(A.A/) 3 + ^(A.A/) 4 + 2Q(X.Aj) S (3.10) 

The stability of the Runge-Kutta schemes is shown graphically in Fig. 3.1 by plotting 
the |G| = 1 contour of Eq. (3.8) for the three-stage, third-order scheme and Eq. (3.10) for 
the five-stage, fourth-order scheme. A selection of XAt in the interior of the closed curve 
yields |G| < 1 , i.e. the scheme is stable. Outside the closed curve, |G| > 1 and the scheme is 
unstable. If the Fourier symbol of the spatial operator, X , is purely imaginary (for example 
the 1-D convection equation) an inspection of Fig. 3.1 reveals that the region, 
-L,<XAt<L,, is stable. If X is purely real (for example the 1-D diffusion equation) the 
region, -L R <XAt<0, \s. stable. The stability limits for these two extreme cases are given in 
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Table 3.3 for the third- and fourth-order Runge-Kutta schemes. The fourth-order scheme 
allows time steps roughly twice that of the third-order scheme. 

A temporal stability analysis of the 3-D convection-diffusion equation is used to 
model the temporal stability of the 3-D Navier-Stokes equations and is postponed until the 
numerical approximations to the spatial derivatives are defined. 


Table 3.3 Stability limits of Runge-Kutta schemes 
for purely imaginary (L/) or real (L R ) spatial operators 


Spatial 

Operator 

third-order, 
three stage 

fourth-order, 
five stage 

Imag, L[ 

1.73 

3.34 

Real, L r 

2.51 

4.65 


3.2 Spatial Discretization 


The numerical approximations to the spatial derivatives appearing in the semi-discrete 
momentum equations, Eq. (3.1), are given in this section. Standard second-order finite dif- 
ference approximations to first derivative suffer from large dispersion errors. Spectral 
methods offer exact differentiation for resolved modes but suffer from high cost and low 
flexibility in that simple domains and boundary conditions are required for their imple- 
mentation. In this study, high-order compact finite differences are used to approximate 
spatial derivatives because of their excellent combination of high-accuracy, flexibility, and 
relatively low cost. 

3.2.1 Numerical Approximation of First Derivative Terms 

The first derivative terms appearing in the governing equations are approximated using 
fourth- and sixth-order compact finite difference schemes proposed by Lele (1992). High 
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accuracy is achieved with relatively small computational stencil sizes by treating the 
derivative terms implicitly: 

where Ax = L X /{N X - \ ),N x '\s the number of grid points, f, represents the first derivative 
of the generic variable <(>, with respect to x, and a , a, b are the coefficients of the compact 
scheme which determine the accuracy. Similar expressions are used for derivatives with 
respect to the y and z directions. For the fourth-order scheme: a = 1/4 , a = 3/2, and 
b - o , and for the sixth-order scheme: a = 1/3 , a = 14/9 , and b = 1/9 . The LHS of Eq. 
(3.1 1) contains the unknown derivatives at grid points i and i+ 1 while the RHS contains 
the known functional values <t>, at the grid points i ± 1 and i ± 2 . 

A comparison of explicit central difference and implicit compact approximations to 
the first derivative is given in Table 3.4. It can be seen that the implicit treatment of the 
derivative results in a smaller or more “compact” stencil for a given order. Also, the lead- 
ing truncation error term for the compact scheme is reduced by 1/4 for the fourth-order 
scheme and 1/9 for the sixth-order scheme compared to explicit central difference 
schemes of the same order. 


Table 3.4 Comparison of explicit central difference 
and implicit compact approximations to the first derivative 


Scheme 

Truncation error 

Stencil Size 

fourth-order central 

ggHfljfflEjjB 

5 

fourth-order compact 


3 

sixth-order central 

(-36/7!)(Ax;)V" 

7 

sixth-order compact 

(-4/7!)(Ax ; )V° 

5 


Writing Eq. (3.1 1) at all grid points results in a tridiagonal system of algebraic equa- 
tions and that is solved efficiently by factoring the LHS into a lower/upper (LU) system 
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once at the beginning of the simulation. The LU factors are stored and then used to solve 
Eq. (3.1 1) for the unknown derivatives. 

The resolution properties *of the numerical approximation to the first derivative are 
analyzed by transforming the 1-D convection equation from physical to wavenumber 
space. In physical space, the 1-D convection equation is given by: 

M-o *«> 

where \u\ is the wave speed which is assumed to be constant. Decomposing into Fourier 
modes as in Eq. (3.6) and evaluating ^ analytically gives: 

|j = (3.13) 

While evaluating the first derivative numerically gives: 


^ <j> (3.14) 

where k* is the numerical wavenumber. For explicit finite difference schemes, a = o, and 
the numerical wavenumber is given by: 


*‘-e2> 


ilkAx 


(3.15) 


l = -N 


while for the tridiagonal compact scheme, the numerical wavenumber is given by: 


* -E 


asin(kAx) + -cos(2A:Ax) 


1 + 2acos(A:Ajc) 


(3.16) 


32 



Note that in general the numerical wavenumber, k*, is complex while the exact wave- 
number, k, is real. For the numerical approximation to yield an exact solution, the follow- 
ing two conditions must be met: 

Real (it*) = k (3.17) 

Imag(£ ) = 0 (3.18) 

It is easy to show that deviations from Eq. (3.17) indicate dispersion errors due to odd 
derivative terms appearing in the truncation error. Deviations from Eq. (3.18) indicate dis- 
sipation errors due to even derivative terms appearing in the truncation error. The real and 
imaginary parts of Eq. (3.16) are plotted separately in Fig. 3.2 for the fourth- and sixth- 
order compact schemes. In addition, some popular explicit numerical approximations to 
the first derivative are plotted for comparison; the standard second-order central difference 
scheme and a third-order upwind biased scheme. 

From Fig. 3.2a, it can be seen that all four approximations do a reasonable job of 
approximating the exact wavenumber (i.e. very low dispersion errors) at very low wave- 
numbers (JfcAx — » 0 ) and that all four approximations do a poor job at very high wavenum- 
bers (JtAjc — » 7 i). For intermediate wavenumbers, the fourth- and sixth-order compact 
schemes provide a much better approximation to the exact wavenumber over a greater 
range of wavenumbers than the explicit schemes. The second-order central difference 
scheme yields a poor approximation to the exact wavenumber for all but the very lowest 
wavenumbers (JtAjt < 0.5 ). From Fig. 3.2b it can be seen that the compact and second-order 
central difference schemes contain no dissipation errors. The third-order upwind scheme 
adds numerical dissipation errors which are largest at high wavenumbers. Spectral meth- 
ods yield exact differentiation for all modes which can be resolved on the specified grid 
and thus correspond to the exact relationship for kAx in Fig. 3.2. 
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Table 3.5 lists some quantitative measures of resolution for the five schemes. The 
wavenumber, k c , defines the region of acceptable accuracy, i.e. for 0 < k < k c , 
k A* — < 0.01 . Modes with k > k c , are not accurately resolved. The quantity, k max , 

defines the maximum value for the modified wavenumber, i.e. for k < k max , the slope of the 
curve is zero. Also listed is the number of spatial grids points per wavelength, 
PPW = 2n/(k c Ax) , to accurately resolve a given mode. From the estimate of PPW, roughly 
five times as many points are required for the second-order central difference scheme to 
achieve the same accuracy as the compact schemes. 


Table 3.5 Resolution measures of various 
numerical approximations to the first derivative 


Spatial Scheme 

Ax 


Points per 
wavelength 

second-order central 

0.22 

1.00 

28.6 

third-order upwind 

0.44 

1.27 

14.3 

fourth-order compact 

1.11 

1.73 

5.6 

sixth-order compact 

1.55 

2.00 

4.1 

spectral 

n 

n 

2 


For non-periodic boundaries, one-sided finite difference expressions are required to 
close the system of equations at the boundary points; i = 1 and i = N for the fourth-order 
scheme and i = 1, 2 and i = N-l, N for the sixth-order scheme. A third-order compact 
boundary scheme is used at i = 1 and i = N with the fourth-order interior scheme: 

3 

(3.19) 

i = 1 

where a bs = 2 and a bSi = -5/2, = 2, a bS} = 1/2 are the coefficients of the third-order 

boundary scheme. A similar equation is used at i = N. 
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For the sixth-order scheme, a boundary and near-boundary scheme are required for 
closure since the interior stencil is pentadiagonal. A fifth-order explicit boundary scheme 
is used at points, i = 1 and i = N: 

4>'i = (3 ' 20) 

i = i 

with coefficients: 


a*,, = -296/105 
a bS2 = 415/48 
a hSi = -125/8 
a bSA = 985/48 


* bSi = - 215/12 
a bj 6 = 791/80 
a bSi = -25/8 
= 145/336 


A different fifth-order explicit near boundary scheme is used at points, i - 2 and i - N-l: 


■>! = < 32l > 
i = 1 

with coefficients: 

a nb) = -3/16 a nbs = 115/144 
a nbi = -211/180 a nbi = -1/3 
a nbi = 109/48 = 23/240 

<> nhi = -35/24 a nfc8 = -1/72 

Similar equations for the boundary and near boundary schemes are used at points i = N 
and i = N-l. 
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3.2.2 Numerical Approximation of the Second Derivative 


The second derivative terms present in the viscous terms of the momentum equation 
and the Laplacian operator of the Poisson equation for pressure are approximated using 
fourth- and sixth-order compact finite differences. Again, higher accuracy is achieved by 
treating the derivative implicitly: 

(Ax) 

+ - 71 (^ 2 -^ + (3.22) 

4(Ajc) 

where <t>", represents the second derivative of the generic variable <t>, with respect to jc, and 
a, a, b are the coefficients of the compact scheme. For the fourth-order scheme: 
a = 1/10, a = 6/5 , and b = 0, and for the sixth-order scheme: a = 2 / 11 , a = 12/11 , and 
b = 3/11 . The tridiagonal system of algebraic equations for the second derivatives are 
solved for in the same manner as the first derivatives. A comparison of explicit central dif- 
ference and implicit compact approximations to the second derivative is given in Table 
3.6. As with the first derivative, the implicit treatment of the second derivative results in a 
smaller stencil size for a given order. The leading truncation error term for the compact 
formulation is reduced by 1/2 for the fourth-order scheme and 1/4 for the sixth-order 
scheme compared to explicit central difference schemes of the same order. 


Table 3.6 Comparison of explicit central difference 
and implicit compact approximations of the second derivative 


Scheme 

Truncation error 

Stencil Size 

fourth-order central 

(-8/6!)(Ax;)> lt ” 

5 

fourth-order compact 

(-3.6/6! )(Ax,)V 6) 

3 

sixth-order cd 


7 

sixth-order compact 

( -(-16.7/8!))(Ax y )V S) 

5 
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For non-periodic boundaries, one-sided finite differences are required to close the sys- 
tem of equations. At i = 1 and i = N, a third-order compact boundary scheme is used: 

, 4 

4>"l + a b$"2 = ~ (3 ' 23) 

,=1 

where a bs = 1 1 and a bSi = 13, a bSi = -27, a bSi = 15, and a bSt = -1 are the coefficients of the 
third-order boundary scheme. For the sixth-order scheme, a near boundary scheme is 
required at i = 2 and i = N-l. The fourth-order interior scheme is used at these points since 
only a three-point stencil is needed. 

A similar analysis of the 1-D diffusion equation is used to investigate the resolution 
qualities of the proposed compact approximation to the second derivative. In physical 
space, the 1-D diffusion equation is written as: 

2 

d<t* 1 d <[> 24) 

Re dx > 

where the term 1/Re represents the diffusion coefficient. Equation (3.24) is transferred to 

2 

wavenumber space by decomposing the solution into Fourier modes. If — \ is evaluated 

dx 

analytically, Eq. (3.24) becomes: 


^ ( 3 . 25 ) 

5 / Re 

Evaluating the second derivative numerically gives: 

| ■ < 3 - 26 , 

where k* is the numerical wavenumber. For explicit finite difference schemes (a = 0), the 
numerical wavenumber is given by: 
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(3.27) 


(*V = 1 




ilk Ax 


The tridiagonal compact scheme represented by Eq. (3.22) yields a numerical wavenum- 
ber: 


(*V = 


(AxY 


2 a[ 1 - cos(/:Ax)] + -[ 1 - cos(2fcAjt)] 


1 + 2acos(A:Ax) 


(3.28) 


Since explicit central and implicit compact difference schemes for the second deriva- 
tive have symmetric stencils ( a, = a_ ( ), the numerical wavenumber is always a real num- 
ber. As a result, there are no dispersion errors in the approximation of the second 
derivative. Only dissipation errors exist. The dissipation errors for the explicit second- 
order central difference, and fourth- and sixth-order compact schemes are shown in Fig. 
3.3 by plotting the numerical wavenumber in Eqs. (3.27) and (3.28). It can be seen that the 
numerical wavenumber for the compact scheme more closely approximates the exact 
wavenumber over a wider range of wavenumbers. Quantitative measures of resolution 
power for the various schemes are given in Table 3.7. It can be seen from the estimate of 
the PPW that roughly twice as many points are required for the explicit second-order cen- 
tral difference to produce the same accuracy as the compact schemes. 


Table 3.7 Resolution measures of various numerical approximations 

to the second derivative 


Spatial Scheme 

(k c Ax) 

(Ca*) 2 

Points per 
wavelength 

second-order central 

0.57 

4.00 

11.0 

fourth-order compact 

1.14 

6.00 

5.5 

sixth-order compact 

1.52 

6.86 

4.1 

spectral 

71 

2 

K 

2 
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Upon discretizing the semi-discrete momentum equation using the compact schemes 


for spatial derivatives outlined above: 


where: 


M + 1 
u i 


u" + b M At[H i -5 x P M ] 


( 3 . 29 ) 


3.2.3 Stability of Runge-Kutta Schemes for the 3-D Convection-Diffusion Equation 

A stability analysis of the 3-D convection-diffusion equation with periodic boundaries 
and uniform grid spacing using the Runge-Kutta scheme is given in this section. The 
results of this analysis are used to model the temporal stability of the 3-D Navier-Stokes 
equations and to select a time step size for the simulations presented in this study. The 
analysis neglects the effects of non-linearity of the convection terms, non-periodic bound- 
aries, the continuity equation, and grid stretching and therefore cannot predict exact stabil- 
ity limits. The approach is to analyze the convection and diffusion equations separately 
and then to combine the two results to determine stability limits for the convection-diffu- 
sion equation. 

Returning to the 1-D convection equation and comparing Eqs. (3.7) and (3.14), it can 
be seen that the Fourier symbol of the spatial operator, is given by, X = -i\u\k . Recall 
that k* is the numerical approximation to the exact wavenumber. The most unstable mode 
in the temporal integration corresponds to the maximum value of XAt over all wavenum- 
bers: 


(kM) max = -i\u\k max At = -iC ^ 


( 3 . 30 ) 
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where C = k max Ax is determined by the numerical approximation to the first derivative and 

is given in Table 3.5 for the various schemes. From the results of the stability analysis of 

the model equation presented in Sec. 3.1, the stability limits for the 1-D convection equa- 

* 

tion are given by: 


-K^iWsn <, 

hi Ax JmaXj 


(3.31) 


where max i denotes the maximum quantity over all grid points. For the 3-D convection 
equation, Eq. (3.31) becomes: 


_i < + 

L, LAx 


M + W* 

Ay Az_L a ,. ; , 


< 1 


(3.32) 


where |v| and |w| are the convection velocities in the y and z directions, respectively. 


Considering the 1 -D diffusion equation, the Fourier symbol of the spatial operator is 
found by comparing Eqs. (3.7) and (3.26). The most unstable mode is given by: 


(*A/) m „ = ~(*‘aJ 2 A/ = -D\— ~ 1 (3.33) 

Re Lj?e(Ax) 2 J 

where D = (k max Ax) 2 is determined by the numerical approximation to the second deriva- 
tive and is given in Table 3.7. The stability limit for the 1-D diffusion equation is given by: 


Dr a t 


< l 


L R L/?e(Ax) 2 

For the 3-D diffusion equation, Eq. (3.34), becomes: 


(3.34) 


DAt 


LpRe 


1 1 1 
— ; + = + ■ 


,1 


< 1 


(3.35) 


Ax)" (Ay)" (A2) 2 -Lx, y , 

The stability limits of the diffusion and convection equation are combined to give the 
stability limits of the convection-diffusion equation. In this case, XAt , possesses both real 
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and imaginary parts. One approach to predicting the stability of the 3-D convection-diffu- 
sion equation, once the computational grid is specified, is to substitute the expressions for 
\At into Eq. (3.8) or Eq. (3.10) and solve the polynomial equation for At at all grid points. 
Selecting the time step as the minimum computed over all the grid points ensures that 
|G| < 1 for all grid points. However, the solution of such a polynomial equation involves 
multiple roots (dependent on the order of the polynomial) and is thus difficult to automate. 

The present approach is to replace the actual contour, Id = 1 , with an ellipse having 
semi-major and -minor axis lengths of Lr and L], respectively. The region inside the 
ellipse is considered to be stable. The stability criterion from this approach is given by: 

A/ < ! (3.36) 

[Litn conv + 

where: 


Lim 


conv 


crM + M + Wl 
LjlAx Ay Az J 


, * - D r i , _l_ . jm 

Lim diff I n Re / A X 2 /a \ 2 /A \ 2 I 

^/? Ae 4Ajt) (Ay) (Az) J 

For free shear flows at moderate to high Reynolds numbers, one expects the stability to 
be governed by the inviscid terms. At lower Reynolds number and/or flows with solid 
boundaries where extremely fine grid spacing must be used, the stability requirements 
could be governed by the viscous terms. An analysis of grids and Reynolds numbers used 
in this study indicates that stability is indeed governed by the inviscid limit and thus an 
explicit time differencing scheme such as the third- or fourth-order Runge-Kutta scheme 
allows reasonably large time steps to be taken. 
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Since there are a number of simplifying assumptions made in extending the stability 
results to the Navier-Stokes equations, the time step used in the simulations is reduced 
from that predicted by Eq. (3.36): 


r r • 2 j , 2, I'* ' ' 

[L/ffl conv + Limjjjf ]nta;t j;t 

where sf is a safety factor, 0 < sf < 1. For the simulations in this study, sf= 0.5 was used. 


In using Eq. (3.36) to estimate the time step requirements for stability, the velocity 
components of the flow field are required to evaluate the frozen convection coefficients, 
M. M, M . In the absence of a representative flow field (such as at the beginning of a simu- 
lation) conservative estimates are used. For non-uniform grids, the local grid spacing 
(Ax ijk , Ay ljk , A z ijk ) is used and it is assumed that the effect of grid stretching on the temporal 
stability is negligible. 


3.3 Enforcement of the Continuity Equation and Poisson Equation for Pressure 


An examination of the governing equations reveals four scalar equations (continuity 
and three scalar components of the momentum equation) in terms of four unknowns (three 
velocity components and pressure). Time derivatives for the velocity components in the 
momentum equation are used to march those equations in time. However, no such time 
derivative exists for pressure, while the continuity equation appears to be an additional 
constraint on the velocity field. The current approach overcomes this problem by taking 
the numerical divergence of the discretized momentum equation and substituting for the 
discretized continuity equation. This results in a Poisson equation for pressure which 
ensures that the velocity field is divergence free at the M+ 1 sub-stage. 
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Applying the divergence operator 5* to the discretized momentum equation gives: 


1 

b M M 


r c / A# + 1 

{M«, 


-«?)} = 5 


( 3 . 38 ) 


The term, 5 x uf ,+ I , represents the discretized continuity equation at the M + / sub-stage 
and is set to zero to enforce the continuity equation. The term, 5 x u * , represents the conti- 
nuity equation at the previous sub-stage M. In practice, this term is retained in the solution 
of the Poisson equation to “kill off’ any accumulating divergence of the velocity field at 
the previous sub-stage. The term, 8 x Hf , is the source term of the Poisson equation and 
represents gradients of the convection and diffusion terms which are known from the pre- 
vious sub-stage. The term, S x 8 X p M , represents the discretized Laplacian operator of the 
pressure. Solving for the Laplacian of the pressure in Eq. (3.38) gives: 


T7 2 M 

V p 


= 8 r 


m n 


Hf + 


b M A/ 


( 3 . 39 ) 


The solution details of the Poisson equation for pressure are given in Chap. 4. 
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Stability diagram for Runge Kutta schemes (dotted line : third-order 
; dashed line : fourth-order). 
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Fig 3.2 Finite differencing error for numerical approximations of the first deriva- 
tive (a) dispersion errors, (b) dissipation errors (zero except for 3rd-order 
upwind). 
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Fig 3.3 



Finite differencing error for numerical approximations of the second 
derivative. 
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Chapter 4 


SOLUTION OF THE POISSON EQUATION FOR PRESSURE 

A significant amount of the total computational time required for the solution of the 
incompressible Navier-Stokes equations (as much as half) is devoted to the enforcement of 
the continuity equation/solution of pressure. This stems from the fact that evolution equa- 
tions exist for the velocity components (i.e. the momentum equations) while none exist for 
the pressure. Instead, an elliptic equation must be solved for the pressure which involves 
the solution of a system of equations and is expensive. Solutions methods for elliptic equa- 
tions generally fall into two categories - direct or iterative. Direct methods usually involve 
some form of Gaussian elimination where the coefficient matrix is first factored into an 
upper and lower matrix and then the solution is computed using back substitution. The 
operation count and memory requirements for this procedure can be prohibitively large for 
the solution of systems involving a large number of unknowns ( - 10 6 in typical 3-D 
problems). The alternative to a direct solution is an iterative procedure where an initial 
approximation to the solution is used to yield an improved solution. This process is 
repeated until the solution is converged within a pre-specified convergence criterion. The 
operation count and memory requirements of most iterative methods are less than that of 
Gaussian elimination. Therefore, the iterative solution procedure is used in this study to 
solve the Poisson equation for pressure. The details of this procedure are outlined in this 
chapter. The performance of the computer code using uniform and curvilinear grids is also 
documented and compared with published computational rates of similar codes. 
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4.1 Discretized Laplacian Operator 


The discrete Poisson equation for pressure which was derived in Sec. 3.3 is given by: 


8 . 8 . 


M 


= 8 . 


H% 


M 1 


b At 


(4.1) 


The LHS of Eq. (4.1) represents a discretized Laplacian operator composed of two 
applications of the first derivative operator, 5 X . It is well known that using two first deriva- 
tive operators to represent the Laplacian operator on non-staggered grids can lead to an 
“odd-even” decoupling of the solution. Indeed, with standard second-order central differ- 
encing for the first derivative operator, the solution at even grid points completely decou- 
ples from the odd grid points, leading to unphysical results. One remedy is to introduce 
terms of the same order as the truncation error which in effect replaces the two first deriv- 
ative operators with a single second derivative operator. This couples the solution at odd 
and even grid points while maintaining the same formal order of accuracy. The Laplacian 
operator is discretized using a single second derivative operator to prevent possible decou- 
pling and Eq. (4.1) becomes: 



= 8 . 


H m + 


M 


b M At 


(4.2) 


where 5 xx p M represents the discrete Laplacian of pressure and is discretized using the 
compact second derivative operator given by Eq. 3.22. 


Writing Eq. (4.2) at all grid points results in a system of equations that is solved at 
each sub-stage of the time advancement scheme. For simplicity, the system of equations 
are defined for the 2-D Poisson equation with periodic boundaries on a uniform grid. The 
Laplacian operator is discretized using the fourth-order tridiagonal scheme defined in Sec. 
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3.2.2. The solution procedure is easily extended to the 3-D Poisson equation. Equation 
(4.2) can be written in the form of a system of equations as: 


AP = [A~ l x B xx + A~yB yy ]P = F 


(4.3) 


A = 


(nj x nj block matrix) 


1 al ol 

al 1 ol 


A yy = 


al 1 al 

[a/ a 7 7 

(nj x nj block matrix) 


Bxx u 1 
h 


(nj x nj block matrix) 


-2/ 7 
7 -27 7 


-21 


B vx = 2 

h y 

1 -2 7 7 

-27 7 -27| 

(nj x nj block matrix) 

p=[p, 

(nj x 1 block vector) 

F =[ F > ’ 

(nj x 1 block vector) 


where l = 


t a 
a 1 a 


a 1 a 
a 1 


(ni x ni scalar matrix) 


where 7 is the ni x ni identity matrix 


R R 


-2 1 1 

n 


1 -2 1 

0 


O 

o 

, where R = 

o 

o 



R 


1 -2 1 

R R 


1 1-2 


(ni x ni scalar matrix) 


where P. = [p, ./> 2 .« 
(ni x 1 scalar vector) 


j 


where V r «ijf 

(ni x 1 scalar vector) 


where P tJ and F Uj = h x {H t + u/(b M At)} t j are the pressure and source term at the ij grid 
point, respectively. The symbols, N x +l,N y +1 denote the number of grid points in the x, y 
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directions, respectively. Values at i = N x +1 and j = N y +1 locations are replaced by values 
at i = 1 andy = 1 since the boundaries are periodic. This change results in a non-zero ele- 
ment in the upper right and lower left comers of the coefficient matrices. The constants, a 
and a, are from the fourth-order tridiagonal scheme for the second deriviative given in Sec. 
3.2.2. The sixth-order compact scheme (with pentadiagonal RHS) can be written in the 
same manner by including the i±2 and j ± 2 terms in the and B yy matrices. 

For non-periodic boundaries, the second derivative boundary scheme given by Eq. 
3.23 is used at the boundary points. In addition, the unknown pressures at the boundaries 
are replaced with their boundary values and those terms are moved and added to the RHS 
of Eq. (4.3). For Dirchlet boundary conditions, such as the freestream conditions given by 
Eq. (2.24), this procedure is straightforward. Neumann boundary conditions, such as those 
applied at the inflow and outflow planes, require that the pressure gradient at the boundary 
be discretized using a first derivative scheme (Eq. 3.19): 



(4.4) 


where the subscripts are used to denote the inflow plane for example. The boundary 
pressure, pjj is then solved for: 



Equation (4.5) is then used to substitute for the boundary pressures in Eq. (4.3). The first 
term on the RHS of Eq. (4.5) is known from the boundary condition and is moved and 
added to the RHS of Eq. (4.3). The second term on the RHS of Eq. (4.5) contains the 
unknown pressures, p 2j and p 3 j , so they are kept on the LHS of Eq. (4.3) and modify the 
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existing terms of the coefficient matrix, A. The resulting system of equations contains only 
the interior pressures as unknowns, Pjj, 2<i<N x -\ and 2<j<N y -\. 

Equation (4.3) results in a “cross” type stencil at the i, j node in which all points along 
lines passing through the central node contribute to the stencil. The coefficients of this 
stencil are implicitly defined in the sense that the matrix operations, [A~ XX B IX + A yy B yy ] , in 
Eq. (4.3) must be performed to determine their values. 

Multiplying Eq. (4.3) by A yy A xx gives: 

[A yy B xx + A yy A xx A-' y B yy )P = A yy A xx F ( 4 . 6 ) 

It is easy to show that the matrices A xx and A yy commute, i.e. A yy A xx = A xx A yy . Using this 
property, Eq. (4.6) simplifies to: 

[A yy B xx+ A xx B yy ]P = A yy A xx F ( 4 . 7 ) 

Using the coefficients of the fourth-order compact approximation to the second derivative 
results in an explicit nine-point, “grid” type stencil for the LHS and RHS of Eq. (4.7): 


p u + 


hv 1 ; 

L W x h]> 

[~ 2 a i % ~ 4 2 )] [/>u+ 1 + P,;_l1 + ' ^)] lP,+ , ’ i ' + + 

+ [P i + i j+ l + Pn.i, j -i + Pi-i,j+] + = F ij + 

L hyZ-i 


( 4 . 8 ) 


a [ F i,j + i +F i,j-l + F i + \,j + F i-\,j] + a l F i+l.j+\ +F i + \, j-\ + F i- 1 , 7+1 +F <- 1 , 7 - 


Using the coefficients of the sixth-order compact approximation to the second deriva- 
tive results in the same nine-point stencil for the RHS and an explicit twenty-one point 
stencil on the LHS: 
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_ 2 r fl+ r \(ux) Pii+ 


*Ahl hV\ 


2tt( b\ a 


2a( b 


+ P U-il + F I-*J * 

y ■* 

[" 1 + P i + 2,j-\ + Pi-2, j+\ + Pi-2, j - \ ] + 

L4 V 

[ 2 I t Pi + \,j + 2 + Pg+ 1,7-2 + 1,7 + 2 + 1 , > - 2 1 + 

+ “ 2 ) f^i+IJ+1 +P i+\,j-\ +P i-l,j+\ + - 1 . > - 1 1 = F ,.j- 


a ^ F i,j* 1 + F i,j- 1 + Fi + \,j + Fi_ i y] + a [F 1 + , y + I + F I+1 y_, + F,_ t J+ 1 + F,_, y _,] 


Thus for uniform cartesian grids, the stencils of Eqs. (4.8) and (4.9) are preferred 
because they require fewer operations and their coefficients are explicit as compared to the 
stencil of Eq. (4.3). The commutative property of the and Ayy matrices is valid even 
with non-periodic boundaries. Numerical experiments confirm that the cross-type stencil 
represented by Eq. (4.3) and the grid-type stencil of Eq. (4.8) or Eq. (4.9) give identical 
results. 


4.2 Iteration Matrix 


A point relaxation scheme is used to iteratively solve the system of equations repre- 
sented by Eq. (4.3) of Eq. (4.8). In this approach, only the value at the central node of the 
stencil, Pjj, is treated as an unknown so that the multi-diagonal system of equations 
degenerates to a diagonal system for one relaxation sweep, which is trivial to solve. This 
process can be written in matrix notation by decomposing the matrix, A, into the sum of 
the diagonal, lower, and upper matrices of A (Briggs 1987): 
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AP = [D-L - U]P = F 


( 4 . 10 ) 


where the matrix, D, is the diagonal matrix of A, and the matrices, L, U, are the lower and 
upper matrices of A, respectively. The solution at the current iteration level, P*> is cor- 
rected with the increment, P\ to yield the solution at the next iteration level, P = P + P' . 
For Jacobi iteration the terms, -LP and -UP, are assumed to be known from the previous 
iterate, P*, and are moved to the RHS so that Eq. (4.10) becomes: 

DP = F +[L+ U]P" ( 4 . 11 ) 

Adding the term -DP* to both sides of Eq. (4.11) casts the equations into incremental 
form: 

DP' = F-AP * = R ( 4 . 12 ) 

where R* is the residual vector. Equation (4.12) defines one iteration of the Poisson equa- 
tion. The increment, P\ is computed by solving the trivial diagonal system of equations 
and is then added to the current iterate to yield the solution at the next iteration level: 

P = p' + D~'r‘ ( 4 . 13 ) 

This procedure is defined as a Jacobi iteration. Weighted Jacobi iteration can under-relax 
or over-relax the iterative process by multiplying the increment, P\ by a relaxation param- 
eter, (o , so that Eq. (4. 13) becomes: 

P = P* + <oD~ l R ( 4 . 14 ) 

where 0<co<i denotes under-relaxation, and w>l denotes over-relaxation. Jacobi 
iteration is equivalent to computing the residual of the current iterate, P*, at all grid points 
followed by an update operation. In this regard, information is held and the solution is 
updated at all grid points simultaneously. Since the computation of the residual vector and 
the updating of the solution vector are completely separate operations, each operation is 
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fully vectorizable. This results in an improved computational rate when those operations 
are performed on vector computers. 

Gauss-Seidel iteration makes one modification to Jacobi iteration - the updated solu- 
tion is used as soon as it is computed. This can be written in matrix form by considering 
Eq. (4.10). The term, -UP, is assumed to be known from the previous iteration and is 
moved to the RHS requiring the solution of the lower triangular system: 

[D-L]P = F+ UP* (4.15) 

In incremental form, Eq. (4.15) becomes: 


P = P* + [D-L)-'r (4.16) 

Moving the term, -UP, to the RHS in Eq. (4. 10) is equivalent to updating the solution vec- 
tor in an ascending order (i.e., Pj j, P 2 j, ■ ■ ■ , P n i,nj )• If the term, -LP, is moved to the RHS 
instead, the variables are updated in descending order (i.e., P n i, n j, P n i.],„j . ■ ■ • , Pi,})- Many 
other updating strategies are possible. For instance, if the grid points are colored black or 
red similar to a checkerboard pattern then during the first sweep all the red points are 
updated, followed by an update of all the black points. This strategy is referred to as red- 
black Gauss-Seidel. The order in which the grid points are visited for Jacobi iteration is 
immaterial since the solution vector is updated only after the residual is computed at every 
grid point. Similar to weighted Jacobi iteration, Gauss-Seidel iteration can be under- or 
over-relaxed by multiplying the increment, P\ by the relaxation factor, to: 

P = P' + w[D - L]~ l R (4.17) 
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4.3 Multigrid Solution 


Relaxation schemes such as Jacobi and SOR applied on a single grid level suffer from 
poor convergence rates as the*number of grid points increase. Through Fourier analysis it 
can be shown (Briggs 1987) that the high frequency error components (i.e. the difference 
between the current iterate and the fully convergenced solution) are removed very quickly, 
while the low frequency components require many iterations to be reduced to an accept- 
able level. In other words, many relaxation schemes are efficient smoothers (i.e. they 
remove high frequency error components with a few iteration sweeps) but are poor solvers 
because they require many iterations to remove low frequency components. 

Multigrid methods overcome these deficiencies by utilizing a hierarchy of grids. 
Smooth error components are transfered to coarser grids where they appear as high fre- 
quency error components and are quickly removed by relaxation sweeps. Relaxation 
sweeps on the coarser grids are also much cheaper to perform. 

A coarse grid correction scheme is utilized in the current study to improve the conver- 
gence rate of the pointwise relaxation scheme on a single grid. Subscripts are used to 
denote grid level, i.e. P h and P 2 h denotes the solution on the fine and coarse grids, respec- 
tively. The symbol, l \ h , is used to denote transfer from the fine to the coarse grid, while 
l h u is used to denote transfer in the opposite direction. The algorithm for one coarse grid 
correction is given below and additional details can be found in Briggs (1987). 

(1) Smooth the current iterate, p ‘ h , on the fine grid v, times: 

A P<>> - F 
A h r h ~ r h 
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(2) Calculate the residual on the fine grid: 


** = Fk-AkK 


,(D 


(3) Transfer (restrict) the residual to the coarse grid, where it is used as the source term 
for the error equation: 




(4) Solve for the error on the coarse grid: 


P'2 h ~ ^2h^2h 

(5) Transfer (prolongate) the error to the fine grid and correct the solution: 

P i 2) = P^ + AhEjh 

(6) Perform v 2 post-relaxation sweeps: 


\P h = P h 

Standard second-order interpolation is used to transfer variables from the coarse to the 
fine grid, while the full weighting operator is used to transfer variables in the opposite 
direction. Although the above algorithm utilizes only two grid levels, improved efficiency 
results from encorporating as many grid levels as possible. In this respect, the direct solu- 
tion of the error equation in step (5) is performed on a very coarse grid requiring a small 
number of operations. Simulations presented in Chaps. 6 and 7 utilize five grid levels. 

Since the simulations are performed on vector computers, Jacobi iteration was utilized 
for relaxation sweeps because it is fully vectorizable. Two pre- and two post-relaxation 
sweeps were performed on each grid. Through numerical experiments, the optimum relax- 
ation factor for the uniform grid formulation was found to be, co = 0.9 . Typically, three 


56 



coarse grid corrections were performed for the solution of the Poisson equation which 
reduced the initial L2 norm of the residual two orders of magnitude. Tests using twice as 
many corrections resulted in a maximum of 3% difference in the pressure. 

4.4 Performance of Computer Code 

The computational rate of the computer code on uniform and curvilinear grids is given 
in this section. The numerical formulation was extended to curvilinear grids which was 
used to concentrate grid points in the mixing layers of the rectangular jet near the inflow. 
Downstream of the potential core, the flow is fully turbulent. Clustering of the grid to 
resolve large gradients due to small scale structures is not possible without the use of time- 
varying, adaptive grids (which is outside the scope of the present study). Therefore, the 
grid is gradually relaxed to uniform spacing in this region. As mentioned in Sec. 4.1, uni- 
form grids enable manipulation of the discrete Poisson equation leading to a reduced sten- 
cil size and cost as compared with the curvilinear grids making the computational rate 
roughly one order of magnitude less. Part of the reason for the increase in computational 
rate is that convergence rates for the multigrid solution of the Poisson equation deteriorate 
with grid clustering and large aspect ratio. Table 4. 1 lists the computational rate for Carte- 
sian and curvilinear formulation. Also listed in the table is the rate of a second-order accu- 
rate formulation from Le and Moin (1994). Therefore, the uniform grid formulation was 
utilized for most of the simulations to be presented in Chaps. 6 and 7. 
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Table 4.1 Computational rate 


Study 

Grid 

Scheme 

Rate 

(cpu/t.s../pt) 

Wilson 

Cartesian 

fourth-order 

compact 

5.7 x 10' 6 

Wilson 

curvilinear 

fourth-order 

compact 

6.2 x 10’ 5 

Le & Moin (1992) 

Cartesian 

second-order 

central 

3.9 x lO' 6 
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Chapter 5 


VALIDATION OF NUMERICAL METHOD 

In this chapter, the numerical formulation is validated through the solution of a wide 
range of benchmark problems. Emphasis is placed on the numerical approximation of spa- 
tial derivatives. In particular, the convection terms (containing first derivatives) present the 
most difficulty in numerical approximation since large dispersion errors exist at high 
wavenumbers ( kAx~n ). It is essential that the numerical approximation to the first deriva- 
tive provide low dispersion errors over a large range of wavenumbers. This is especially 
true in 3-D simulations where reducing the required number of grid points by half in each 
coordinate direction leads to eight times fewer total grid points. 

In Chap. 3, the theory showed that compact schemes require roughly five times fewer 
points to accurately resolve a given mode compared to the standard second-order central 
difference approximation to the first derivative. This theory is tested by solving some prac- 
tical problems ranging from the 1-D convection equation to the 2-D Navier-Stokes equa- 
tions. 

The purpose of the present chapter is to (i) demonstrate the resolution qualities of the 
compact schemes and (ii) compare and contrast results from the compact scheme with 
results from other popular schemes. 
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5.1 1-D Convection Equation 


The first problem to be solved is the 1-D convection equation, Eq. (3.12), which tests 
the time advancement scheme and the numerical approximation to the first derivative. The 
exact solution corresponds to convection of the scalar profile at the constant wave speed, c. 
Distortion in the shape of the profile indicates dissipation and/or dispersion errors in the 
solution. The convection of a Gaussian profile was solved using three approximations to 
the first derivative; (i) a second-order central difference, (ii) a third-order upwind, and (Hi) 
the fourth-order compact approximation outlined in Section 3.2. 1 . The third-order Runge- 
Kutta scheme was used to advance the equation in time for all spatial schemes. In addition, 
the CFL number (and thus the time step) was kept small so that resulting errors are due to 
the spatial formulation. 

The parameters and initial conditions are those proposed at the ICASE/LARC Work- 
shop on Benchmark Problems in Computational Aeroacoustics, Hardin et al. (1995): 


u(x, 0) = 0.5exp 



(5.1) 


-20 < x < 450 , N x = 470 grid points, c = 1 

Since the specified grid is relatively coarse in comparison with the initial conditions, this 
problem provides an excellent test of the resolution power of the numerical approxima- 
tion. Figure 5.1 shows the computed solutions at / = 400 after the profile has convected to 
x = 400. There is little discernible difference between the exact solution and the solution 
with the fourth-order compact scheme. However, the solutions with the second-order cen- 
tral difference and the third-order upwind approximations show greatly reduced peak val- 
ues and large, dispersive waves trailing the Gaussian profile. The errors from the second- 
order central difference scheme are the most severe. 
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It is difficult to determine by inspection what portion of the error is dispersive and 
what portion is dissipative. The solutions are transformed into wavenumber space using a 
Fourier transform method and compared with the exact solution in Fig. 5.2 to address this 
issue. The graph displays the resulting complex Fourier coefficient in polar form with the 
amplitude displayed in Fig. 5.2a and the phase angle in Fig. 5.2b. It can be seen from Fig. 
5.2a that the solutions computed with the second-order central difference and fourth-order 
compact schemes predict the correct amplitude for all modes. The amplitude of the solu- 
tion computed with the third-order upwind scheme is reduced or dissipated, especially at 
higher wavenumbers. Figure 5.2b shows that the fourth-order compact scheme predicts 
the correct phase angle even for the highest wavenumbers. 

The phase angle from the second- and third-order solutions are only correctly pre- 
dicted for the very lowest wavenumbers (k < 0.2 for the second-order solution and k < 0.3 
for the third-order solution). Large dispersion errors are evident at high wavenumbers. The 
above trends in the numerical solutions are consistent with the dissipation/dispersion error 
theory for the 1-D convection equation and show the resolution power of the compact 
schemes. 

The second problem, also proposed at the ICASE/LARC Workshop on Benchmark 
Problems, is the solution of the 1-D convection equation in a spherical coordinate system. 
The governing equation takes the form: 


^ + “+^ = 0 (5.2) 

dt x ox 

5 < x < 450 , N x = 445 grid points 


Initial Conditions: 


«(*, 0 ) = 0 
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Boundary Conditions: 


u(5, t) = sin^-'j 

Figure 5.3 shows the exact solution at t = 400 which corresponds to a damped sine wave 
due to the addition of the u/x term in the governing equation. Fig. 5.4 shows computational 
results for the region, 200 < x < 220 using the third-order upwind approximation to the 
first derivative on two grids and the fourth-order compact scheme. The solution with 1 6 
points per wavelength shows a greatly reduced amplitude and a phase shift relative to the 
exact solution. It takes roughly 64 PPW (not shown) to reproduce the exact solution with 
the third-order upwind approximation. The fourth-order compact approximation is able to 
reproduce the exact solution with 8 PPW. 

5.2 2-D Convection Equation 

Multidimensional effects of the numerical formulation are explored by solving for the 
convection of an inverted cone around a circle. This problem is governed 2-D convection 
equation: 


du du du _ 
di +Cx fa + C >ty ~ ° 


(5.3) 


where c x = -y and c y = x, are the convection speeds in the x and y directions, respectively. 
The initial conditions are that of an inverted sharp cone centered at x, y = -0.5, 0. The exact 
solution corresponds to the cone being convected counterclockwise in a circular path of 
radius, r Q = 0.5 with a period of 2 % . Distortion of the shape of the cone is an indication of 
dispersion and/or dissipation errors. 


Figure 5.5 shows computed results after one revolution of the cone using ( a ) a third- 
order upwind approximation and (b) a fourth-order compact approximation to the first 
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derivatives on a 32 x 32 grid with uniform spacing. This grid defines the shape of the cone 

with a maximum of 8 points in each coordinate direction. The exact shape of the cone is 

included to the right of the computed solution at x,y = 0.5, 0 for comparison purposes. The 

% 

third-order solution (Fig 5.5a) shows that the sharp point of the cone is greatly diffused 
and that dispersion errors are evident trailing the cone. A grid of 128 x 128 (or 32 points 
defining the shape of the cone) must be used with the third-order upwind approximation 
before the shape of the cone is faithfully reproduced. The fourth-order compact solution 
(Fig 5.5b) shows that the shape of the cone is not distorted as it is convected around the 
circle on the 32 x 32 grid. Indeed, the only noticeable error is a very small “grid to grid” 
oscillation due to the absence of physical viscosity in this problem and numerical viscosity 
in the compact scheme. 

Figure 5.6 shows results for the same problem after one revolution obtained by Orszag 
(1971) using (a) second-order Arakawa finite-difference, (b) fourth-order Arakawa finite- 
difference, and (c) a spectral approximation to the first derivatives on a 32 x 32 grid. The 
finite difference solutions show errors similar to the third-order solutions in Fig. 5.5. The 
spectral method, which provides exact differentiation for all wavenumbers representable 
on the 32 x 32 grid, convects the cone without distorting its shape. Thus, the solution using 
the compact scheme is closer to the spectral solution than the solutions obtained with con- 
ventional finite difference schemes. The higher accuracy and resolution characteristics are 
achieved by the implicit treatment of the derivative. Even though the stencil size of the 
compact scheme is finite, the implicit treatment of the derivatives makes the scheme glo- 
bal much like that of spectral methods. 
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5.3 2-D Euler/Navier-Stokes Equations 


In the previous sections, the effect of numerical approximation on the accuracy of the 
convection terms was documented. In this section, the accuracy of the enforcement of the 
continuity equation through the solution of the Poisson equation for pressure is 
documented by solving the 2-D Euler/Navier-Stokes equations. Since the Navier-Stokes 
equations contain viscous terms, the numerical approximation to the second derivative is 
also tested. The test problems chosen for validation contain many features of the 3-D jets 
which are simulated in the current study. In this respect, the test problems are not merely 
academic exercises. There benchmark problems are solved; (i) a temporally-developing 
plane mixing layer (2-D Stuart’s problem), (ii) 2-D viscous wave decay, and (iii) the 
doubly periodic jet. Problems (i) and (ii) have exact solutions while results of (iii) will be 
compared to a highly resolved spectral simulation. 

5.3.1 2-D Temporally-Developing Plane Mixing Layer 

Exact solutions to the Euler or Navier-Stokes equations for general flows do not exist 
due to the non-linearity of the convection terms. However, under special conditions exact 
solutions may be found. An exact solution for the temporally-developing mixing layer was 
first published by Stuart (1967). The initial conditions for the 2-D Stuart’s problem corre- 
spond to a steady hyperbolic tangent function for the streamwise velocity component with 
a periodic array of vortex cores in the mixing region which cause the solution to vary in 
time. The wavelength of the disturbance corresponds to the neutral mode such that the dis- 
turbance is convected in the streamwise direction with no change in amplitude. The exact 
solution for the streamwise velocity component, u, and the transverse velocity component, 
v is given by: 


u(x, y,r) = c + 


Csinh(y) 

Ccosh(y) + Acos(jr - ct ) 
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( 5 . 4 ) 


_ Asin(jc-cf) 

v{x,y,t) - (j cos j,(j,) + Acos(x-ct) 

where A = Jc 2 - 1 is a parameter which controls the strength of the perturbation and c is 

the convective speed of the mixing layer. The flow is periodic in the streamwise direction 

with length, L x = 2n , 0 < x < 2n . The flow is infinite in the transverse direction but in this 

study is truncated at a finite distance, -L y <y< L y , such that the zero-traction freestream 

boundary condition outlined in Section 2.2.3 is well approximated. Tests which vary the 

transverse domain height, 2 Ly, show that Ly — 10 is sufficiently large to implement this 

boundary condition. The exact solution is shown in Fig. 5.7a with parameters, c = 1 , A = 1/ 

2. A uniform, cartesian grid is used for the simulations in this section. Unless otherwise 

specified, the third-order Runge-Kutta scheme is used for time advancement and time 

steps are sufficiently small so that spatial errors are dominant. 

Figure 5.7b shows the numerical solution at t = 207c (ten flow through times) on a rel- 
atively coarse grid of 13 (streamwise) x 41 (transverse) using the fourth-order compact 
approximation of convection terms and pressure. The solution of pressure involves the 
computation of the source term and the discretization of the Laplacian operator in Eq. 
(3.39). In addition, once the Poisson equation is solved for pressure, the gradient of pres- 
sure is computed which is required to advance the momentum equation in time. Therefore, 
the phrase “fourth-order solution of pressure” corresponds to the source and pressure gra- 
dient terms computed with the compact first derivative scheme outlined in Section 3.2.1 
while the Laplacian operator is discretized using the compact second derivative scheme 
outlined in Section 3.2.2. 

Even though the grid is relatively coarse (13 streamwise points per wavelength and 
roughly 8 points in the mixing region at y~0), there is little discernible difference 
between the exact and numerical solutions after ten flow through times. It is important to 
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check the convergence of the error as the grid is refined to expose any coding errors, to 
demonstrate that the order of error convergence seen in practical computations is that pre- 
dicted by a Taylor series analysis, and to gain confidence in the numerical formulation. 
Tables 5. 1 and 5.2 give a quantitative measure of the L2 and maximum errors in the veloc- 
ity components at t = 0.1 using the fourth- and sixth-order compact approximations to the 
convection terms and solution of pressure, respectively. Solution errors from three grids 
are shown where the grid spacing in the x and y directions is halved from coarsest to finest 
grid. The results show that the L2 and maximum errors converge at roughly the rate pre- 
dicted by a Taylor series analysis as the grid is refined. The order, N, is computed using the 
solution error from three grids of spacing, h, 2h, and 4h: 



where <|> A , <j> 2/1 , <J» 4A are the errors on the h, 2/i, and 4h grids, respectively. In using Eq. (5.5) it 
is assumed that the solution is fully resolved on all three grids and that the leading trunca- 
tion error term is dominant (Demuren and Wilson 1994). 

Table 5.1 Solution errors for 2-D Stuart’s Problem at / = 0.1 using fourth-order 
compact approximation for convection terms and solution of pressure. 


Grid (ni x ni) 

Max U Error 

Max V Error 

L2 Norm U 

L2 Norm V 

13x41 

■jSEDESH 

0.21 x 10' 2 

0.18 x 10' 3 

■EBDESH 

25x81 



besbemh 

ESH 

49 x 161 

HsUQisyH 

Msmm 


kszqeh 

Order (N) 

4.2 

4.1 

4.5 

■ 


Table 5.2 Solution errors for 2-D Stuart’s Problem at t — 0.1 using sixth-order 
compact approximation for convection terms and solution of pressure. 


Grid (ni x ni) 

Max U Error 

Max V Error 

L2 Norm U 

L2 Norm V 

13x41 

Hi&DQfll 


0 97 X 10‘ 4 

WESESEtBI 

25x81 

oa 


bzesemb 

■ 

49 x 161 

HZE2XESH 

0.18 x 10' 5 



Order (N) 

5.4 

5.7 

6.3 

6.2 
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To address the effect of computing the pressure with a lower-order formulation, the 2- 
D Stuart’s problem was solved using second-order central, fourth-order compact and 
sixth-order compact approximation of the convection terms but a second-order central dif- 
ference solution of the pressure. The results of the three computations are shown in Tables 
5.3 - 5.5. The results of the three computations show that the lower-order solution of pres- 
sure results in the overall convergence of the error being second-order, even if the convec- 
tion terms receive a higher-order treatment. All terms must be discretized using higher- 
order approximations to achieve higher-order error convergence rates. 

Table 5.3 Fourth-order compact approximation for convection terms/ 
second-order solution of pressure. 


Grid (ni x ni) 

Max U Error 

Max V Error 

L2 Norm U 

L2 Norm V 

13x41 

IKESDSsiHi 

0.11 x 10' 1 

0.20 x 10 ^ 

0.20 x 10^ 

25x81 

HE11QQH! 

IQgQQBI 

0.42 x 10' 5 

0.44 x I0* 3 

49 x 161 


HiHIlQH 

o.n x io'* 

0.11 x 10‘ 3 

Order (N) 

2.2 

2.2 

2.4 

2.2 


Table 5.4 Sixth-order compact approximation for convection terms/ 
second-order solution of pressure. 


Grid (ni x ni) 

Max U Error 

Max V Error 

L2 Norm U 

L2 Norm V 

13x41 

^gjgQQjgB 

0.11 x 10' 1 

0.20 x 10 ^ 

0.20 x lO*- 4 

25 x 81 

HEISSI2H 




49 x 161 


mmmm 

HumiEyH 

mSBQHl 

Order (N) 

2.1 

22 

2.3 

2.2 


Table 5.5 Second-order central difference approximation for convection terms/ 

second order solution of pressure. 


Grid (ni x ni) 

Max U Error 

Max V Error 

L2 Norm U 

L2 Norm V 

13x41 


0.15 xlO" 1 

0.21 x 10 -!J 

HiEiHEMH 

25x81 

HsUIIIsSjfll 



rnia 

49 x 161 



0.11 x I0' 3 

0.13 x 10' 3 

Order (N) 

2.4 

1.9 

2.3 

2.1 


The solution of the 2-D Stuart’s problem validates the numerical formulation for the 
enforcement of the continuity equation and the solution of the Poisson equation for pres- 
sure. In addition, it has been shown that the zero-traction freestream boundary condition 
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for shear flows is a valid approximation provided that the freestream boundary is located a 
sufficiently far distance from the mixing region. 

5.3.2 Viscous Wave Decay . 

The numerical treatment of viscous terms is validated by solving the 2-D viscous wave 
decay problem which is governed by the Navier-Stokes equations. The domain for this 
problem is periodic in both the x and y directions where periodic boundary conditions are 
applied. The exact solution is given by: 

-(-) 

u(x, y, t) = -cos(jt)sin(y)e ^ Re ' 

-f-1 

v(x, y, t) = sin(jc)cos(y)e 

where Re = 20, L x = Ly = 1 . The exact solution consists of sinusoidal waves in the x and y 
directions which decay in time. Table 5.6 shows the L2 norm of the error at t = 0.025 using 
the fourth- and sixth-order compact approximations for convection and diffusion terms 
and the solution of pressure. The results are compared to the fourth-order, Essentially 
Non-Oscillatory (ENO) scheme from Weinan and Shu (1992). The error converges at 
fourth- and sixth-order rates thus validating the numerical treatment of the viscous terms 
and again validating the convection terms and the solution of pressure. The error of the 
ENO scheme converges at a fourth-order rate, but is more than two orders of magnitude 
greater than the fourth-order compact results. The error magnitude of the sixth-order com- 
pact formulation on the 128 x 128 grid has reached the round-off error level ( - 10"' 3 ) of 
the Cray supercomputer, indicating that extremely accurate results are obtained on aver- 
age-sized grids. 
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Table 5.6 Solution errors for 2-D viscous wave decay. 


Grid (ni x nj) 

4th oa compact 

6th oa compact 

3rd(4th) oa ENO 

16 x 16 


0.10 x 10 7 

- 

32x32 

0.77 x 10‘ 8 

0.15 x 10’ 9 

0.53 x 1 0' 5 

64 x 64 

0.47 x 10' 9 

0.27 x 10 " 

0.32 x 10‘ 6 

128 x 128 

0.71 x 10 10 

0.11 x 10‘ 12 

0.20 x 10’ 7 

Order(N) 

4.0 

6.0 

4.0 


5.3.3. Doubly Periodic Jet 

The last validation test is the solution of the doubly periodic jet which is governed by 
the 2-D Euler equations. The initial conditions correspond to a jet or “top hat” profile for 
the streamwise velocity component. The initial conditions and problem parameters are: 


u(x, y, 0) 


r, anh p-s2)] ;yS „-| 


v(jc, y, 0) = 8sin(x) 

where L x = L y = 2n, 6 = 0.05, and p = n/15. The flow begins with two parallel, finite- 
thickness shear layers, one with positive vorticity and the other with negative vorticity. A 
small amplitude perturbation is provided through the transverse velocity component which 
causes the shear layers to roll up into vortex cores as they evolve. Between the vortex 
cores, the shear layers are stretched and thinned as they are wound around the vortex 
cores. Without viscosity, the scales of motion become smaller and smaller so that eventu- 
ally resolution is lost on any fixed grid. This problem represents a worst case scenario 
where the scales of motion cannot be resolved on the grid. Therefore, it is an extremely 
demanding test of the resolution power of the numerical formulation. 


69 





















Figure 5.8 shows vorticity contours of the evolution of the doubly periodic jet using 
spectral methods on a highly resolved grid of 512 x 512 from Weinan and Shu (1992). An 
18th-order filter has been used to remove energy at the highest wavenumbers which would 
otherwise contain aliasing errors. Since spectral methods yield exact differentiation and no 
unphysical oscillations occur during the time interval, 0 < t < 8, we can infer that this 
simulation is extremely accurate during this period. The vorticity contours at t = 10 show 
the beginning of unphysical oscillations or “wriggles” which is a sign that the 
computations are under-resolved and the smallest scales of motion are determined by the 
grid and not the physics of the problem. 

Figure 5.9 shows vorticity contours for an unfiltered simulation of the doubly periodic 
jet using the sixth-order compact formulation for convection terms and the solution of 
pressure on a 128 x 128 grid. Even though the simulation contains 1/1 6th the number of 
grid points used in the spectral simulation, all of the relevant features of the jet are 
captured. The results show that the location of the braids and details of the vortex cores are 
well predicted. However, for t> 6, the results are not as smooth as the spectral simulation 
and unphysical “wriggles” appear in the braid region at t = 10. 

Figure 5.10 shows vorticity contours from a simulation using the sixth-order compact 
scheme where the velocity field is filtered every ten time steps with an explicit sixth-order 
filter. For this computation, the velocity field is filtered in the x and y coordinate directions 
to remove energy at the highest wavenumbers ( kAx-n ). The unfiltered quantity, 0, is fil- 
tered in the x coordinate direction to produce the filtered quantity, <j> : 


i + <{>, + a<(> 1+ 1 (5.6) 

= a<t>, + 1(4>, + , + i ) + 2 + - 2 ) + \ (♦,- + 3 + <t>, - 3) 
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where a, a, b, c, and d are the coefficients of the filter. The coefficients for the sixth-order 
explicit filter are summarized in Table 5.7 along with a second-order explicit and fourth- 
order compact filter. A similar equation exists for filtering in other coordinate directions. 
Since the boundaries for this problem are periodic, no additional boundary schemes are 
required. For problems with non-periodic boundaries, such as the simulations presented in 
Chap. 7, the following fourth-order boundary filters are used at the first three grid points: 

$1 = n<t>i + 16 (4<,)2 ~ 6<t, j + 4<t>4 - ♦s) 

3 1 

4*2 = ^<t>2 + 15^1 + 6<t>3 -4<t>4 + <t> 5 ) 

h = + 4«t> 2 + 4*4 - <(. 5 ) (5.7) 

By transforming the filter from physical to wavenumber space, its effect on the various 
modes can be clearly shown. The transfer function defines the filtering operation in wave- 
number space (Lele 1992): 


T(kAx) - 


a + bcos(kAx) + ccos(2A:Ax) + dcos(3kAx) 
1 + 2acos(fcAx) 


(5.8) 


The transfer function for the sixth-order explicit filter is plotted in Fig. 5. 1 1 along with 


a second-order explicit and fourth-order compact filter for comparison. The LES presented 
in Chap. 7 utilizes the fourth-order compact scheme together with this fourth-order 
compact filter. It can be seen that the fourth- and sixth-order filters eliminate the highest 
wavenumber ( it Ax - n ) while leaving the low wavenumbers unchanged. The fourth-order 
compact filter is a high-bypass filter in the sense that relatively high wavenumber modes 
(it Ax < 1.5) pass through the filter without being changed. Recall from Fig. 3.2 that the 
approximation to the first derivative is accurate for the low to moderately high 
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wavenumbers and that wavenumbers, kAx~n, are not well represented. Therefore, the 
effect of the filter is to remove energy from those modes which are not well represented. 

The vorticity contours from the filtered simulation (Fig. 5.10) show that the filtering 
operation removes the unphysical oscillations while maintaining the fine scale details of 
the flow. In Fig. 5.12, a simulation using the fourth-order ENO scheme from Weinan and 
Shu (1992) is shown for comparison purposes.The simulation using the ENO scheme with 
the same grid shows that the braids and vortex cores are diffused and that information is 
lost for t > 8 . A simulation of the doubly periodic jet at t = 10 using the sixth-order com- 
pact scheme and explicit filter is shown in Fig. 5.13 on a 256 x 256 grid. 


Table 5.7 Coefficients for filters 


Scheme 

a 

a 

b 

C 

d 

second-order, explicit 

0 

1/2 

1/2 

0 

0 

fourth-order, compact 

0.475 

(5 + 6a) 
8 

0 +2a) 
2 

(l-2a) 

8 

0 

sixth-order, explicit 

0 

11/16 

15/32 

-3/16 

1/32 
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Fig 5.1 Solution to the 1-D convection equation at t = 400 in physical space for 
various finite difference approximations of the first derivative term. 


73 




0.0100 



wavenumber, k 



wavenumber, k 


Fig 5.2 Solution to the 1-D convection equation at / = 400 in wavenumber space 
for various finite difference approximations of the first derivative term, 

(a) amplitude and (b) phase angle. 
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Fig 5.5 Numerical solution of the rotating cone problem after one revolution on a 
32 x 32 grid (a) third-order upwind scheme, (b) fourth-order compact 
scheme. Numerical solution is shown to the left, exact solution to the right. 
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Fig 5.10 
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Numerical solution of the double shear layer using sixth-order compact 
scheme with filtering on a 128 2 grid, (a) t = 4, (b) t = 6, (c) t = 8, and 
(d) t = 10. 
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Fig 5.11 Transfer function for filter schemes. 




Fig 5.12 Numerical solution of the double shear layer using fourth-order ENO 
scheme on a 128 2 grid from Weinan and Shu (1992), (a) t = 4,(b)t = 6, 

(c) t = 8, and (d) t = 10. 
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Chapter 6 


DIRECT NUMERICAL SIMULATION - RESULTS 

Results from the direct numerical simulation of rectangular jets at Reynolds number of 
Re De = 750 are presented in this section. The Reynolds number of the jet flow is necessar- 
ily low so that all the scales of motion can be resolved with reasonable computational 
resources and no subgrid scale model is needed. The effect of Reynolds number on jet 
dynamics is addressed in the next chapter where large-eddy simulations are preformed at 
Re De = 75,000. Emphasis is placed on LES of non-circular jets (Chap. 7) at higher Rey- 
nolds numbers where the flow is turbulent as occurs in experiment and nature. Spatial sim- 
ulations are performed in this study where the domain is fixed in space and the flow enters 
and exits the domain through inflow and outflow boundaries, respectively. This is in con- 
trast to a temporal simulation where the domain moves in the streamwise direction such 
that a small region of the flow is followed in time. The spatial reference frame is the one 
that occurs in nature and is thus preferable. It is also much more computationally demand- 
ing because the entire domain of interest must solved be simulated for all times. 

6.1 Discrete Mode Forcing 

As discussed in Sec. 2.2.1, time dependent boundary conditions are applied at the 
inflow to model the jet nozzle a short distance downstream of the exit ( x/D e < 1) and to 
promote the development of coherent structures within the computational domain. In this 
section, results are presented for simulations using discrete modes given by Eq. (2.19). 
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Simulations with discrete modes are used to study shear layer instabilities and the effects 
of forcing. Two cases are simulated with discrete modes corresponding to; (i) the funda- 
mental mode, (to = 0.22) at an rms intensity of 3% of the mean core velocity, U a , and (ii) 
the fundamental mode (to = 0.22) and the first subharmonic mode (to = 0.11) at 1.5% 
intensity each. Cases (iii) and (iv) are presented in the next section and utilize the broad 
mode forcing at 3% and 15% total rms intensity, respectively. The domain and grid dimen- 
sions for the four runs are summarized in Table 6.1. 


Table 6.1 Summary of parameters for DNS 


Case 

Forcing function 

L x x LyX L z 

N x xN y xN z 

(i) 

fund. @ 3% 

1 

I 

| 

(H) 

fund. + 1st Subh. 
@ 1 .5% each 

12 x 10 2 

129 x 129 2 

(Hi) 

broad @3% 

12 x 10 2 

129 x 129 2 

(iv) 

broad @15% 

12 x 10 2 

129 x 129 2 


The domain length, L x , listed in Table 6.1 includes a buffer layer of two diameters in 
length, giving an active computational domain of ten diameters. The streamwise grid spac- 
ing results in 10 points per fundamental wavelength. In Chap. 4 it was shown through the 
solution of several benchmark problems, that 8 points per wavelength provides adequate 
resolution with the fourth-order compact scheme. 

Contours of vorticity magnitude for case (i) at t = 2 flow through times (the time to 
travel from inflow to outflow at the convective velocity) are shown in Figs. 6. la-e. Con- 
tours in the major and minor axis planes (Figs. 6.1a and b ) reveal that the shear layers near 
the inflow plane roll up at the fundamental frequency within the first diameter. For the 
region, x/D e > 5, the fundamental mode has saturated and decayed. With no additional har- 
monic modes present in the forcing function, the jet does not transition to turbulence at 
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this Reynolds number. Near the outflow plane, x/D e ~ 10, the jet width in the minor axis 
(x-y) plane becomes larger than that in the major axis (x-z) plane indicating that axis 
switching has taken place. Figs. 6.1 c - e shows the vorticity magnitude at three cross-sec- 
tional planes ( x/D e = 0, 5, and 10), while Figs. 6.1/- h shows streamwise vorticity at three 
cross-sectional planes. ( x/D e = 2.5, 5, and 10). Even though there is no streamwise vortic- 
ity introduced at the inflow plane (not shown), by x/D e = 2.5 (Fig. 6.1 f) four pairs of 
streamwise vortices have formed in the higher curvature comer regions. This results in the 
distortion of the initially rectangular boundary layer at x/D e = 0 to the diamond pattern at 
x/D e = 10. 

Vorticity contours for case (ii) alt = 2 flow through times show that the effect of add- 
ing the first subharmonic mode to the fundamental mode is to promote the formation of 
additional structures for the region, x/D e > 5. In addition, the shear layer is dominated by 
intense vortices in the minor axis plane spaced at the first subharmonic wavelength. By 
examining the cross-sectional profiles in Fig. 6.2c -e it is apparent that rapid shear layer 
growth takes place in the major axis plane and that axis-switching does not take place. 
Cross-sectional contours of streamwise vorticity at x/D e = 5 (Fig. 6.2g) show that four pair 
of streamwise vortices develop very close to the jet centerline which leads to a a partial 
bifurcation of the jet near x/D e = 10. The streamwise vortices must be generated from the 
redistribution of azimuthal vorticity since they are not present at the inflow. 

6.2 Broad Mode Forcing 

Broad mode forcing is utilized to model naturally developing (unforced) jets with tur- 
bulent boundary layers. The resulting forcing function is somewhat random and does not 
contain symmetries present in discrete mode forcing. Figure 6.3 shows representative vor- 
ticity magnitude contours alt = 2 flow through times. 
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In contrast with simulations performed with discrete mode forcing, the shear layers do 

not roll up in a periodic fashion. Instead, non-symmetric, random structures are formed 

which are characteristic of naturally-developing jets. Contours at the outflow (Fig. 6.3f) 

♦ 

show that the jet width in the minor axis plane is larger than that in the major axis plane, 
indicating that axis-switching takes place near the outflow plane. 

The impact of the unsteady structures on the mean flow is examined by time-averaging 
the results of case (iii) with broad mode forcing at 3% intensity over a period of 8 flow 
through times. The results from the first flow through time are not included in the time- 
averaging so that transients resulting from the initial conditions are convected out of the 
domain. After the transient period, the flow field at every other grid point and every fourth 
time step is saved to disk. The results are then post-processed to compute statistical quan- 
tities such as first and second moment quantities, two-point correlations, and budgets. This 
procedure results in roughly 1 188 samples in the period of 8 flow through times. The ade- 
quacy of sample size is addressed in the next chapter. 

The time-averaged jet major and minor axis widths, entrainment, decay of centerline 
velocity, and fluctuating centerline velocity from case (iii) are shown in Fig. 6.4 as a func- 
tion of streamwise coordinate. The results at this level of intensity show that significant 
unsteadiness does not occur until, x/D e - 7 . This results in a small growth of the jet widths 
and no distinctive end to the potential core region as observed in higher Reynolds number 
experimental jets at x/D e = 4-5. Decay of centerline velocity plots from the DNS of a rect- 
angular jet at Repg = 800 of Miller et al. (1995) also reveal no distinctive end of the poten- 
tial core within their computational domain of nine diameters. 

In order to test the effect of forcing intensity, a fourth simulation is presented at a 
higher intensity of 15%. Jet widths in the major and minor axis planes, entrainment rate. 
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centerline velocity, and fluctuating velocity are given in Fig 6.5. Detailed profiles of time- 
averaged velocity, pressure, and Reynolds stress components for case ( iv) are shown in 
Figs. 6.6 - 6.15. Time-averaged and fluctuating velocity components are normalized with 
the local velocity on the jet centerline, Uql and transverse coordinates are normalized by 
the local jet width defined by the transverse distance where U/Ucl = 0.5. The location, Y 0 , 
Z 0 = 0, 0 denotes the position of the jet centerline. 

Figure 6.5d shows that the effect of increasing the intensity of the inflow forcing func- 
tion to 15% is to promote unsteadiness at x/D e ~ 5 . The centerline velocity also begins to 
decay slowly at roughly x/D t - 5 . Figure 6.6 shows time-averaged streamwise velocity 
along the minor and major axis. The experimental results of a jet issuing from a contoured 
rectangular nozzle of AR = 2 at Re De = 10 5 (Quinn 1995) are also shown. This experimen- 
tal jet shows axis-switching at x/D r - 12 . The profiles from the DNS at x/D e =9 show good 
agreement with experimental profiles in the minor axis plane only. Fluctuating velocity 
profiles shown in Figs. 6.10-6.12 are generally underpredicted in comparison with the 
experimental profiles at higher Reynolds number. This is not surprising since the center- 
line fluctuating velocity is still increasing near the end of the computational domain of ten 
diameters (Fig. 6.5d) while the experimental value reaches its peak at x/D e = 5-6 and then 
levels off. 


89 




Fig 6.1 Contours of vorticity magnitude (a) - (e) and streamwise vorticity (f) - (h) 
for case (i) at t = 2 flow through times for fundamental forcing function, (a) 
minor axis plane, z/D e = 0, (b) major axis plane, y/D e = 0, (c) cross flow 
plane, x/D e = 0, (d) x/D e = 5, (e) x/D e = 10, (f) x/D e = 2.5, (g) x/D e = 5, and 

(h)x/D e = lO. 
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Fig 6.2 Contours of vorticity magnitude (a) - (e) and streamwise vorticity (f) - (h) 
for case (ii) at t = 2 flow through times for fundamental and first subhar- 
monic forcing function, (a) minor axis plane, z/D e = 0, (b) major axis plane, 
y/D e = 0 ,(c) cross flow plane, x/D e = 0, (d) xlD e - 5, (e) x/D e = 10, (f) x/D e = 
2.5, (g) x/D e = 5, and (h) x/D e = 10. 
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10 D„ 



Fig 6.3 Contours of vorticity magnitude for case (iii) at t = 2 flow through times 
for broad mode forcing function, (a) minor axis plane, zJD e = 0 ,(b) major 
axis plane, y/D e = 0, (c) cross flow plane, x/D e =2.5, (d) x/D e = 5, (e) x/D e = 
7.5, and (f) x/D e = 10. 
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Jet Width 




Fig 6.4 Time-averaged quantities versus streamwise distance for DNS of rectan- 
gular jet with broad mode forcing at 3%; (a) jet widths (solid - major axis 
plane, dashed - minor axis plane), (b) entrainment ratio, (c) decay of cen- 
terline velocity, (d) fluctuating velocity (solid - u rm JU q, dashed - u rm J 

U C L )• 
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Fig 6.5 Time-averaged quantities versus streamwise distance for DNS of rectan- 
gular jet with broad mode forcing at 15%; (a) jet widths (solid - major 
axis plane, dashed - minor axis plane), (b) entrainment ratio, (c) decay of 
centerline velocity, (d) fluctuating velocity (solid - u rms /Uo, dashed - u rms / 

Uci)- 
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(Z-Z 0 )/Z , /2 


Time-averaged streamwise velocity, U, for the DNS of rectangular jet, (a) 
minor axis and (b) major axis (solid - experiment of Quinn 1995, x/D = 10). 
x/D = ; + : 0.0 ; *: 1.5 ; Q : 3.0 ; O • 4 - 5 5 A : 6 -° ; □ : 7 * 5 5 X : 9 -°- 










(Z-Z 0 )/Z 1/2 


Time-averaged lateral velocity, V, for the DNS of rectangular jet, (a) minor 
axis and (b) major axis. For symbols see Fig 6.6. 



















(Z — Z 0 )/2 1/2 

Fig 6.10 Time-averaged fluctuating streamwise velocity, u rms , for the DNS of rect- 
angular jet, (a) minor axis and (b) major axis (solid - experiment of Quinn 
1995, x ID = 10). For symbols see Fig 6.6. 








(Z-Z 0 )/Z , /2 


Fig 6.11 Time-averaged fluctuating lateral velocity, v rms , for the DNS of rectangu- 
lar jet, (a) minor axis and (b) major axis (solid - experiment of Quinn 
1995, x/D = 10). For symbols see Fig 6.6. 







(Z-Z 0 )/Z 1/2 


Fig 6.12 Time-averaged fluctuating transverse velocity, for the DNS of rect- 
angular jet, (a) minor axis and (b) major axis (solid - experiment of Quinn 
1995, x/D = 10). For symbols see Fig 6.6. 















(Z-Z 0 )/Z , /2 


Fig 6.15 Time-averaged Reynolds shear stress, <vV’>, for the DNS of rectangular 
jet, (a) minor axis and (b) major axis. For symbols see Fig 6.6. 
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Chapter 7 


LARGE EDDY SIMULATION-RESULTS 

This chapter discusses the results of large eddy simulation (LES) of non-circular jets at 
a higher Reynolds number, Re De = 75,000. At this Reynolds number, it is not practical 
with current computer resources to fully resolve all scales of motion of the jet flow. The 
Smagorinsky model outlined in Sec. 2.1.2 is included in the computation to account for 
unresolved scales. In addition, the velocity field is filtered every five time steps using the 
fourth-order compact filter outlined in Sect. 5.3.3. Filtering removes the highest wave- 
number mode which is not accurately resolved, while passing through modes of low to 
moderate wavenumber. Results are presented for jets with initially rectangular and ellipti- 
cal cross-sections with low aspect ratio, AR = 2. In addition, LES of a circular jet is per- 
formed to quantify the effects of non-uniform azimuthal curvature present in the non- 
circular jets. 


7.1 Rectangular Jet 


7.1.1 Simulation Parameters 

In Sec. 2.2.1, the mean velocity of the rectangular jet at the inflow plane was specified. 
Broad mode forcing functions are used to promote unsteadiness and to model the naturally 
developing jet observed experimentally. A low level of forcing is used at this Reynolds 
number which corresponds to a maximum intensity at the inflow plane of u rm /&x = 0 . 03 . 
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The domain dimensions are 12 a: 10 at 10 and grid dimensions are 129 3 . The last two 
diameters of streamwise domain length are used as a buffer domain giving an “active” 
computational length of 10 diameters in the streamwise directions. The grid is distributed 
uniformly in all three directions. 

7.1.2 Effect of Grid Resolution 

The effect of grid resolution on the LES results is presented in this section. In the LES 
approach, it is often difficult to separate the effect of truncation errors, aliasing errors, and 
SGS model errors as the grid is refined. This is contrasted with the DNS approach where 
all scales are resolved such that no SGS modeling errors are committed and the truncation 
error reduces at a known rate with grid refinement. Therefore, it is easier to isolate the 
effects of the numerical formation from those due to the physics of the problem. 

The approach taken in the current study is to solve the problem numerically on a very 
coarse grid (denoted as the level 1 grid) using the same parameters identified in the previ- 
ous sections. After roughly ten flow through times, the results on the coarsest grid are pro- 
longated to a grid with double the number of grid points in each direction. Those results 
are used as initial conditions for a simulation on the level 2 grid which is run for ten flow 
through times. This process is repeated until it is no longer practical to double the number 
of grid points due to computer resource limitations. This procedure has the advantage of 
efficient use of computer resources in that solutions on the coarsest grid are very cheap 
due to the small number of grid points and large time step taken. Possible problems with 
the processing of the runs and the numerical formulation are identified with minimal com- 
puter resources, and confidence is gained in the solution procedure by the time the finest 
(and most expensive) grid level is reached. The prolongated result from the coarser grid 
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level also provides a “better” initial condition than starting from an artificially created ini- 
tial condition which is inconsistent with the governing equations. 

A by-product of this procedure is that results are generated on a number of grids with 
increasing resolution which can be used to address the issue of the sensitivity of the results 
to grid resolution. A summary for the grids used in the resolution study is given in Table 
7.1. along with the number of streamwise points per fundamental wavelength. The 
increased resolution for the first three grid levels is obtained by doubling the number of 
grid points while keeping the domain dimensions constant. Continuing this trend to the 
fourth grid level would require 208 x 256 2 grid points which is very computationally 
expensive. A sufficiently long simulation on this grid would require roughly 350 Cray C- 
90 hours and 920 megawords of memory. An alternative is to examine the resolution of the 
first five diameters (on the fourth grid level only) by fixing the number of grid points and 
halving the domain dimensions. This region includes the entire potential core and the thin 
shear layer region near the inflow where resolution requirements are greatest. 


Table 7.1 Summary of grids used in resolution study 


Grid Level 

Ax, Ay, Az 

N r N y N z 

Lx, Ly L z 

Points per wave 

1 




2.5 

2 

0.193 x0.156 2 

52 x 64 2 

10 3 

5 

3 

0.096 x 0.078 2 

104 x 128 2 

10 3 

10 

4 

0.048 x 0.039 2 

104 x 128 2 

5 3 

20 


Figure 7.1 shows the jet width, decay of centerline velocity, and fluctuating centerline 
velocity for the 4 grid resolutions. The results show that with the level 1 grid with roughly 
2.5 points per fundamental wavelength, very little unsteady is resolved by the grid (see 
Fig. 7.1c). This results in a small spread rate in the minor axis plane and no axis-switching 
(lower curves, Fig. 7.1a). Increasing the grid resolution to the second level increased the 
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unsteadiness to 5% at x/De = 10. The results on the third level show a clear end to the 
potential core, followed by a leveling off of the unsteadiness, axis-switching is clearly 
shown in Fig. 7.1a. Increasing the grid resolution to the fourth level has a small impact on 
those trends. The biggest difference between the level 3 and 4 results is that the end of the 
potential core is shifted upstream by roughly 0.7 diameters (from 4. ID on level 3 to 3.4D 
on level 4). To show this point, the results on the third and fourth grid levels are replotted 
in Figs. l.\d - /with the x/De coordinate shifted 0.7 diameters downstream for the fourth 
level. The graph shows that the results on the third level are relatively insensitive to grid 
refinement after the correction for potential core length. The discrepancy in the potential 
core length is most likely due to differences in the broad mode forcing function on the 
third and fourth grid level. The maximum rms intensity of the streamwise velocity pertur- 
bation is fixed at 3% of the core velocity for all grids. However, since the time step is 
halved with each grid refinement, the spectra content of the broad mode forcing is neces- 
sarily different on each grid. As a result of the grid resolution study, simulations in Chaps. 
6 and 7 utilize the grid resolution of grid level 3. Future work should include a simulation 
on the 10 3 domain using 208 x 256 2 grid points, although it will be costly. 

7.1.3 Velocity Spectra 

Figure 7.2 displays time traces and power spectra for u, v, w, and p at six spatial loca- 
tions throughout the computational domain. The first three locations correspond to the 
inflow plane ( x/De = 0), roughly halfway (x/De = 4.13) and near the outflow plane (x/De = 
8.34). The transverse location for the first three points is fixed in the center of the major 
axis shear layer (y, z = 0.31, 0), while the last three transverse locations are fixed in the 
minor axis shear layer, (y, z = 0, 0.63). The power spectra is computed by transforming the 
time traces from the temporal domain to the frequency domain, and are useful in determin- 
ing the range of temporal scales and dominant frequencies. The results in Figs. 7.2 a and 


108 



1.2d show that the forcing function is composed of a wide band of frequencies of similar 
strength. At five diameters (Figs. 1 .2b and e) the power spectra show that a dominate fre- 
quency has emerged, while the highest frequency modes are several orders of magnitude 
smaller. Near the outflow (Figs. 7.1c and f) the range of scales is more broad suggesting 
that the flow is turbulent. A line with slope,/' 5/3 , is included in the figures which is used to 
infer the presence of an inertial subrange and turbulent flow. The power spectra near the 
outflow (Figs. 7.2c and f) show roughly one decade with the/ scaling. 

7.1.4 Instantaneous Flow Field 

Contours of streamwise velocity, pressure, vorticity magnitude, and streamwise vortic- 
ity at t = 9 flow through times are shown in Figs. 7.3 - 7.9. Negative contours are drawn 
with dashed lines while positive contours are drawn with solid lines. Local maximum and 
minimum values are also indicated in the figure. Contours in the minor axis plane (Fig. 
7.3c) show that the shear layers roll up at roughly, x/De = 3, similar to a plane mixing 
layer. Pressure contours (Fig. 7.3b) show that the shear layers above and below the jet cen- 
terline roll up in an organized and staggered fashion. At x/De = 4-5, the unsteady struc- 
tures from the upper and lower shear layers meet at the jet centerline (y/D = 0), thus 
signaling the end of the potential core. Downstream, the flow is characterized by smaller 
scale, less organized structures. For x/De < 4 the vorticity is dominated by the azimuthal 
components, i.e. there is no streamwise vorticity present. However, for x/De > 4, it is obvi- 
ous that significant streamwise vorticity has been generated. 

Figure 7.4c shows that shear layer roll up in the major axis plane is suppressed by roll 
up in the minor axis plane. Thus, the spreading of the jet in the major axis plane is sup- 
pressed resulting in a switching of the major and minor dimensions of the jet at x/De = 7. 
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The axis-switching at t = 9 is clearly seen by comparing streamwise velocity contours at 
the cross-section, x/De = 2.44 (Fig. 7.6) and x/De = 9.75 (Fig. 7.9). 

7.1.5 Time-Averaged Flow Field 

In the previous section, the flow field was examined at one instant in time to describe 
the formation and development of the unsteady structures. In this section, the impact of the 
unsteady structures on the mean flow is examined by time-averaging the results over a 
period of 1 1 flow through times. The results from the first flow through time are not 
included in the time-averaging so that transients resulting from the initial conditions are 
convected out of the domain. The results are then post-processed to compute statistical 
quantities. This procedure results in roughly 1634 samples in the period of 1 1 flow 
through times. The adequacy of sample size is addressed in the next section. The raw data 
files from the simulation occupy roughly 1 6 gigabytes of disk space. 

Contours of Time-Averaged Flow Field 

The time-averaged contours of streamwise velocity are shown in Fig. 7. 10. Comparing 
contours of streamwise velocity confirms that axis-switching has indeed taken place at 
roughly seven diameters as was suggested by examining the instantaneous contours at t = 
9 flow through times. Figure 7.10e and f shows that roll up and interaction of structures 
result in rapid spreading in the minor axis plane only. The end of the potential core occurs 
at roughly x/De = 4.5, where the velocity along the jet centerline is no longer equal to the 
core velocity and begins to decrease due to entrainment of ambient fluid. 

Profiles of First and Second Moment Quantities 

In this section detailed profiles of time-averaged velocity, pressure and the six unique 
Reynolds stress components are shown. In addition, jet widths in the major and minor axis 
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planes, entrainment rate, centerline velocity and fluctuating centerline velocity are given. 
The numerical results are compared with available experimental results. 

Time-averaged and fluctuating velocity components are normalized with the local 
velocity on the jet centerline, U CL and transverse coordinates are normalized by the local 
jet width defined by the transverse distance where U/Dq^ = 0.5. Figure 7.11a shows that 
axis-switching takes place at x/De = 6.8. Tsuchiya (1985) reported the axis-switching 
location of their 2: 1 AR jet from a smoothly contoured nozzle to be x/De = 6.0. A second 
axis-switching was reported at x/De = 25 which is beyond the computational domain used 
in this study. Experiments of turbulent jets from contoured nozzles (Quinn 1995) show the 
first axis-switching to be as far as x/De = 12. A measure of entrainment into the jet (Fig. 
7.11b) is provided by computing the difference in mass flux at the cross-section, x y and 
that at the inflow, Q(x) - Q 0 , where Q(x) = JJp U(x,y,z)dydz and Q 0 = jjpU(0,y,z)dydz. 
Figure 7.1 lc shows the decay of centerline velocity with streamwise distance. The end of 
the potential core is predicted at x/De = 4.0, which is in excellent agreement with the 
experimentally measured value of x/De = 4.1 by Tsuchiya (1985). The development of 
fluctuating velocity on the jet centerline is shown in Fig. 7.1 Id. For x/De < 4 the velocity 
fluctuations on the centerline are small. Near the end of the potential core, the fluctuating 
velocity rises sharply which is a result of the unsteady structures from the mixing layers 
meeting at the centerline. Downstream of the potential core, the fluctuating velocity is 
roughly constant when scaled with the local centerline velocity. 

Profiles of time-averaged streamwise velocity are shown in Fig. 7.12 at various 
streamwise locations. For comparison, the experimental results of Quinn (1995) at x/De = 
10 are included. The LES results at x/De = 9 are in excellent agreement with the experi- 
mental results. Self-similar profiles are predicted in the minor axis plane downstream of 
the potential core, while non-similar ones are predicted in the major axis plane. This 
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observation is consistent with the experimental studies of Sforza et al. (1966), Trentacoste 
and Sforza (1967), and Tsuchiya (1985). Figure 7.13 shows some self-similarity of the 
time-averaged lateral velocity component, V, along the y axis at the last two cross-sections 
( x/De > 7). Note that negative (positive) transverse velocity at the top (bottom) edge of the 
jet is consistent with fluid being entrained into jet. The lateral velocity along the z axis 
should be zero due to symmetry. The numerical results predict small lateral velocity along 
the z axis. Time-averaged transverse velocity profiles are shown in Fig. 7.14 which are 
clearly non-similar along the z axis. The transverse component should be close to zero 
along the y axis due to symmetry. Krothapalli et al. (1981) reported similar trends for rect- 
angular jets from smoothly contoured nozzles at higher aspect ratios AR = 16.7. Time- 
averaged pressure profiles displayed in Fig. 7.15 show small positive values in the poten- 
tial core followed by large negative values downstream of the potential core. 

The normal components of the Reynolds stress tensor (plotted as rms values) are 
shown in Figs. 7.16 - 7.18, while the shear components are shown in Figs. 7.19 - 7.21. In 
general, profiles of the normal components are relatively flat along the y axis near the cen- 
terline, while off center peaks are present along the z axis downstream of the potential 
core. Upstream of the end of the potential core, the peaks in the normal components corre- 
spond to unsteadiness in the mixing layers separating the jet core from the ambient sur- 
roundings. 

Peak values are higher (20 - 30%) than those reported for a rectangular jet by Quinn 
(1995). Similar trends exist for the shear components. The shear components of the Rey- 
nolds stresses show that the dominant components are predicted for <u’v’> along the y 
axis (Fig. 7.19a) and <u ’w’> along the z axis (Fig. 7.20b). 


112 



7.1.6 Two-Point Velocity Correlations 


In this section, results of two-point velocity correlations taken on the jet centerline at 
two cross-sections are presented. The correlations are used to estimate the size of the 
coherent structures along the jet centerline downstream of the potential core. In addition, 
the correlations are used to access the placement of the freestream boundaries in the com- 
putation. Large correlation coefficients at the edges of the computational domain indicate 
that freestream boundaries are too close to the centerline. 

The transverse size of the coherent structures was estimated using two-point spatial 
correlations which requires prescribing the spatial separation, r, along the transverse direc- 
tion with zero time separation. The two-point correlation is given by: 


R.jiK + » = 


t) 2 )(uj{> 0 + >, r) 2 )] 1/2 


( 7 . 1 ) 


where < > denotes time-averaging, > n - (x n , y a , z„) denotes the location of the correlation, 
and f denotes the separation distance and direction of the correlation. An estimate of the 
structure size can be obtained by fitting a parabola through the data points near y = 0 and z 
= 0 in Figs. 7.22 - 7.24. In each figure, the top plot represents the correlation with the sep- 
aration distance along the y axis, while the bottom plot is along the z axis. The correlations 
show that events near the jet centerline are uncorrelated with those near the edge of the 
domain, thus justifying the placement of the free stream boundaries in the computation. 


7.1.7 Budgets 


In this section a detailed term by term budget is presented for the resolved time-aver- 
aged momentum equations and the Reynolds stress transport equations. The imbalance in 
the terms of the resulting equation can then be used to access the adequacy of the sample 
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size. The momentum and Reynolds stress budgets can also be used for turbulence model- 
ing of complex three-dimensional jets. Computation of budget terms for the plane mixing 
layer was performed by Rogers and Moser (1994). Demuren et al. (1996) used DNS data 
to develop turbulence models for the pressure diffusion in mixing layers and wakes. 

Momentum Equation 

The mean momentum equation for the resolved velocities is derived by time-averaging 
the filtered equations of motion given by Eq. (2.8). Omitting the details of the derivation, 
the result becomes: 


du, rr du, BP 1 d z U i \ , „ (B\ x 

Bt ~ U Q Xj Bx, Redxjdxj 


( 7 . 2 ) 


where (/,■ = («,) and P = (<j>) denotes the time-averaged velocity and pressure, respec- 
tively. The first three terms on the RHS of Eq. (7.2) represent the convection, pressure gra- 
dient, and viscous diffusion terms, respectively. 


For statistically stationary results, the terms on the RHS of Eq. (7.2) should sum to 
zero indicating that the time derivative of the average velocity is also zero. The fourth term 
on the RHS of (7.2) represents the resolved stresses due to unsteadiness of the velocity 
field and is similar to the Reynolds stress term in the Reynolds Averaged Navier-Stokes 
(RANS) equations. This term is computed directly in the large-eddy simulation approach. 
The last term on the RHS of Eq. (7.2) represents the contribution of the unresolved 
stresses. The concept of LES is that a large portion of the energy containing large scale 
stresses (fourth term on RHS) is resolved in the computation, compared with the unre- 
solved portion which is modeled (last term on RHS). 
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Profiles of budget terms for the U momentum equation are shown in Fig. 7.25 along 
the y and z axis at the location, x/De = 9. Since the unresolved stresses, q^, are not avail- 
able from the LES results, the eddy viscosity model given by Eq. (2.9) is used to compute 
the last term in Eq. (7.2). The resulting profiles in Figs. 7.25 - 7.30 have been passed 
through a filter to remove some high frequency noise which tends to obscure the trends. 
The sum of the terms on the RHS of Eq. (7.2) is represented by the symbols. The profiles 
show that the imbalance is relatively small compared to the convection and resolved Rey- 
nolds stress terms, and that the convection terms balance the Reynolds stress terms. The 
results also show that the unresolved SGS stresses are quite small compared with the 
resolved Reynolds stresses, validating the LES. The budgets for the V momentum equa- 
tion in Fig. 7.26 show that along the y axis, the pressure gradient term balances the 
resolved Reynolds stress term at y/D - ±2 . 

The budget terms along the z axis are smaller in comparison as would be expected 
from symmetry arguments. Figure 7.27 shows similar trends for the W equation. The sense 
of the pressure gradient term for the V and W equation is consistent with fluid being 
entrained into the jet on a time-averaged basis. 

Resolved Reynolds Stress Equation 

The transport equations for the resolved Reynolds stresses are derived by first subtract- 
ing the filtered equations of motion, Eq. (2.8), from the time-averaged filtered equations, 
Eq. (7.2). The result is a transport equation for the i th component of the resolved fluctuat- 
ing velocity component. Time-averaging the quantity; uj ■ [u/equation] + u- ■ [uj equation] 
gives the transport equation for the resolved Reynolds stresses. Again, omitting the details 
of the derivation: 
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( 7 . 3 ) 



d , , , 
-d7 k (Ui u * u * 



■/£♦<■&> 


2 ^duj | i 3 2 <«,-'«/) 
Re'dx k dx k Re dxjdxj 



The first six terms on the RHS of Eq. (7.3) represent convection of the Reynolds 
stresses by the mean flow, turbulent production, turbulent dissipation, fluctuating velocity/ 
pressure gradient coupling, turbulent diffusion, and viscous diffusion, respectively. The 
last term on the RHS of Eq. (7.3) represents the contribution of the unresolved stresses to 
the resolved Reynolds Stresses equation. As with the momentum budget, the SGS model 
is used to estimate this term. 


Figure 7.33 shows profiles of the <u’u’> budget of the resolved Reynolds stress equa- 
tion at the same location as the momentum budgets. It can be seen that a large positive pro- 
duction term is opposed by a negative velocity pressure gradient term. Convection and 
turbulent diffusion are also significant. Turbulent dissipation and viscous diffusion are 
small. The contribution of the SGS term is estimated to have a modest negative contribu- 
tion. Budgets for the <u’u’> and <u’u’> are shown in Figs. 7.29 and 7.30a, respectively. 

The results show that the imbalance is of the same order as some of the individual 
terms near the jet centerline. The imbalance is most likely due to the unresolved terms 
which are estimated. The dominate role of the very smallest scales (which are not resolved 
in LES) is in the dissipation of turbulent kinetic energy into heat. Since the smallest scales 
are not resolved, the turbulent dissipation term in Eq. (7.3) does not resemble the term in 
the fully resolved Reynolds stress equation. As a result, there is a large positive imbalance 
due to this discrepancy. Indeed, in some experimental studies which measure budgets of 


116 



the Reynolds stress equations, the imbalance is labeled as the turbulent dissipation 
because it cannot be accurately measured. The budget for the <u’u’> equation from the 
DNS of the rectangular jet (case iii of Chap. 6) is shown in Fig. 7.30b. The overall balance 
of terms is better since the dissipation scales are resolved. 

7.2 Elliptic Jet 

In this section, the results of a large eddy simulation of a 2:1 AR elliptic jet at Re De = 
75,000 with broad mode forcing are presented. The grid and domain parameters are iden- 
tical to those described for the LES of the rectangular jet. The jet boundary layer at the 
inflow is generated by using, n = 2 for the exponent of the super-elliptic coordinate sys- 
tem. Figure 7.31 shows that the mixing layer in the potential core of the elliptic jet rolls up 
preferentially in the minor axis plane at the expense of mixing in the major axis plane. The 
roll ups are not as orderly and well-defined as in the rectangular jet. Contours at the cross- 
section, x/D e = 4.88 (Fig. 7.35d) reveal the generation of significant streamwise vorticity 
around the circumference of the jet boundary layer resulting in the distortion of the elliptic 
cross-section. It is apparent that the major and minor axis have already switched by the 
cross-section, x/D e = 9.75 (Fig. 7.37). 

Time-averaged contours of streamwise velocity shown in Fig. 7.38 reveal that axis- 
switching also takes place in the elliptic jet. Time-averaged results for the elliptic jet simu- 
lation are presented in Figs 7.39 - 7.49. The instantaneous results are averaged over eight 
flow through times resulting in a sample size of roughly 1200. The plot of jet widths (Fig. 
7.39a) confirms that axis-switching has taken place at x/D e = 5.9 and that the width in the 
major axis plane actually decreases slightly in agreement with the experimental results of 
unexcited 2:1 AR elliptic jets studied by Hussain and Husain (1989). The axis-switching 
location in that study occurred at roughly x/D e = 5.0. Fig. 7.39b shows that the elliptic jet 
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entrains more fluid than the rectangular jet which is consistent with experimental observa- 
tions. A summary of axis switching location is provided in Table 7.2 along with those 
from experiment. 


Table 7.2 Summary of axis switching location from LES 


Study 

Geometry 

AR 

Location, x/D e 

Tsuchiya et al. (1985) 

Rectangular 

2 

6.0 

Quinn (1985) 

44 

44 

11.0 

Current 

44 

44 

6.8 

Hussain and Husain (1989) 

Elliptic 

44 

5.0 

Current 

(4 

44 

5.9 


7.3 Round Jet 


Large eddy simulation of a round jet is performed to provide a control for the simula- 
tions of rectangular and elliptic jets. Parameters for the round jet simulation are identical 
to those for the non-circular jets, with the obvious geometrical difference. For the round 
jet simulation, some asymmetry of profiles along the y and z axis exist beyond 6 diameters 
due to insufficient sample size. Perfectly symmetrical time-averaged profiles are difficult 
to achieve in experiments and computations of fully turbulent flow. 

Contours of instantaneous streamwise velocity, pressure and vorticity for the round jet 
are shown in Figs. 7.50 and 7.51 in the x-y and x-z plane, respectively. Cross-sectional 
contours are shown in Figs. 7.52 - 7.56. The contours in the x-y and x-z plane show that 
regular, planar rings are not formed with the broad mode forcing function. Cross-sectional 
contours at x/D e - 4.88 (Fig. 7.54) show significant generation of streamwise vorticity 
along the jet boundary layer resulting in a distortion of the initially circular profile. Further 
downstream, the vortices breakdown into small scale structures. 
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Time-averaged contours of streamwise velocity shown in Fig. 7.57 reveal that jet 
widths in the x-y and x-z planes grow at the same rate and axis-switching does not take 
place. Time-averaged quantities from the round jet simulation are shown in Fig. 7.58 for 
comparison with the non-circular jet simulations. 
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Fig 7.1 Effect of grid resolution, (a) -(c) time-averaged quantities versus x, (d) - (f) 
level 4 corrected for length of potential core; (dot-dash : level 1 grid ; 

3 dots-dash : level 2 grid ; solid : level 3 ; dash : level 4). 
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Fig 7.2b Time traces (first column) and power spectra (second column) for U, V, W, 
and P on minor axis, (x, y, z) = (4.13, 0.31, 0) Re De = 75,000. 














































Fig 7.3 Contours in the minor axis plane {zJD = 0) at t = 9 for LES of rectangular 
jet; (a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.4 Contours in the major axis plane (z/D = 0) at t = 9 for LES of rectangular 
jet; (a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.5 Contours at the cross-section, x/D = 0, at t = 9 for LES of rectangular jet; 

(a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.6 Contours at the cross-section, x/D = 2.44, at t = 9 for LES of rectangular 
jet; (a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.7 Contours at the cross-section, x/D = 4.88, at t = 9 for LES of rectangular 
jet; (a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.8 Contours at the cross-section, x/D = 7.31, at t = 9 for LES of rectangular 
jet; (a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.9 Contours at the cross-section, x/Z) = 9.75, at t = 9 for LES of rectangular 
jet; (a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.11 Time-averaged quantities versus streamwise distance for LES of rectangu- 
lar jet; (a) jet widths (solid - major axis plane, dashed - minor axis plane), 
(b) entrainment ratio, (c) decay of centerline velocity, (d) fluctuating veloc- 
ity (solid - u rms /U 0 , dashed - u rms /U CL ). 
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Fig 7.12 Time-averaged streamwise velocity, U, for the LES of rectangular jet, (a) 
minor axis and (b) major axis (solid - experiment of Quinn 1995, x/D = 10). 
For symbols see Fig. 6.6. 


136 




(Z-Z 0 )/Z 1/2 


Fig 7.13 Time-averaged lateral velocity, V, for the LES of rectangular jet, (a) minor 
axis and (b) major axis. For symbols see Fig 6.6. 
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Fig 7.14 Time-averaged transverse velocity, W, for the LES of rectangular jet, (a) 
minor axis and (b) major axis. For symbols see Fig 6.6. 
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Fig 7.16 Time-averaged fluctuating streamwise velocity, u rms , for the LES of rect- 
angular jet, (a) minor axis and (b) major axis (solid - experiment of Quinn 
1995, x/D = 10). For symbols see Fig 6.6. 
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Fig 7.17 Time-averaged fluctuating lateral velocity, v rms , for the LES of rectangular 
jet, (a) minor axis and (b) major axis (solid - experiment of Quinn 1995, x/ 
D - 10). For symbols see Fig 6.6. 
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Fig 7.18 Time-averaged fluctuating transverse velocity, for the LES of rectan- 
gular jet, (a) minor axis and (b) major axis (solid - experiment of Quinn 
1995, x/D = 10). For symbols see Fig 6.6. 
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Fig 7.19 Time-averaged Reynolds shear stress, <u’v’>, for the LES of rectangular 
jet, (a) minor axis and (b) major axis. For symbols see Fig 6.6. 
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Fig 7.25 Budget terms for the U -momentum equation at x/D = 9.38 for LES of 
rectangular jet, (a) minor axis and (b) major axis. 
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Fig 7.27 Budget terms for the W -momentum equation at x ID - 9.38, for LES of 
rectangular jet (a) minor axis and (b) major axis. 
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Fig 7.29 Budget terms for the <v V> equation at x/D = 9.38, for LES of rectangular 
jet (a) minor axis and (b) major axis. 









Fig 7.30b Budget terms for the <u , u , > equation at x ZD = 9.38, for DNS of 

rectangular jet (case iii, Chap. 6) (a) minor axis and (b) major axis. 
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Fig 7.31 Contours in the minor axis plane ( z/D = 0) at t = 3 for LES of elliptic jet; 

(a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.32 Contours in the major axis plane {z/D = 0) at t = 3 for LES of elliptic jet; 

(a) streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) 
streamwise vorticity. 
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Fig 7.33 Contours at the cross-section, xJD = 0, at t = 3 for LES of elliptic jet; (a) 
streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.34 Contours at the cross-section, yJD = 2.44, at t = 3 for LES of elliptic jet; (a ) 
streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.35 Contours at the cross-section, x/D = 4.88, at t = 3 for LES of elliptic jet; (a) 
stream wise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.36 Contours at the cross-section, x/D 
stream wise velocity, (b) pressure, 
wise vorticity. 
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Fig 7.37 Contours at the cross-section, x/D = 9.75, at t = 3 for LES of elliptic jet; (a) 
streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.39 Time-averaged quantities versus streamwise distance for LES of elliptic 
jet; (a) jet widths (solid - major axis plane, dashed - minor axis plane), (b) 
entrainment ratio, (c) decay of centerline velocity, (d) fluctuating velocity 
(solid - u^Uq, dashed - u rms /U CL ). 
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Fig 7.40 
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Time-averaged streamwise velocity, U, for the LES of elliptic jet, (a) minor 
axis and (b) major axis (solid - experiment of Quinn 1995, x/D = 10). 
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Fig 7.42 Time-averaged transverse velocity, W, for the LES of elliptic jet, (a ) minor 
axis and (b) major axis. For symbols see Fig 6.6. 






Fig 7.43 Time-averaged pn 
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Fig 7.44 Time-averaged fluctuating stream wise velocity, u ms , for the LES of elliptic 
jet, (a) minor axis and (b) major axis (solid - experiment of Quinn 1995, x! 
D = 10). For symbols see Fig 6.6. 
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Fig 7.45 Time-averaged fluctuating lateral velocity, v rms , for the LES of elliptic jet, 
(a) minor axis and (b) major axis (solid - experiment of Quinn 1995, x/D = 
10). For symbols see Fig 6.6. 
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Fig 7.47 Time-averaged Reynolds shear stress, <uV>, for the LES of elliptic jet, 
(a) minor axis and (b) major axis. For symbols see Fig 6.6. 
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Fig 7.48 Time-averaged Reynolds shear stress, <u’ w’>, for the LES of elliptic jet, 
(a) minor axis and (b) major axis. For symbols see Fig 6.6. 
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Fig 7.49 Time-averaged Reynolds shear stress, <vV>, for the LES of elliptic jet, 
(a) minor axis and (b) major axis. For symbols see Fig 6.6. 
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Fig 7.50 Contours in the X-Y plane (z/D = 0) at / = 3 for LES of round jet; (a) 
streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.52 Contours at the cross-section, x/Z) = 0, at t = 3 for LES of round jet; (a) 
streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.53 Contours at the cross-section, x/D = 2.44, at/ = 3 for LES of round jet; (a) 
streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.54 Contours at the cross-section, x/D = 4.88, at t = 3 for LES of round jet; (a) 

streamwise velocity, (b) pressure, (c) vorticity magnitude, and (d) stream- 
wise vorticity. 
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Fig 7.58 Time-averaged quantities versus stream wise distance for LES of round 

jet; (a) jet widths (solid - major axis plane, dashed - minor axis plane), 
(b) entrainment ratio, (c) decay of centerline velocity, (d) fluctuating 
velocity (solid - u^/Uq, dashed - u rms /U CL ). 
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Chapter 8 


SUMMARY AND CONCLUSIONS 


Three-dimensional simulations of turbulent jets with rectangular and elliptical cross- 
section were simulated with a newly developed numerical formulation. At low Reynolds 
numbers the full Navier-Stokes equations were solved, while at higher Reynolds numbers 
the filtered equations of motion were solved along with a sub-grid scale model. The time- 
dependent results from the simulation are used to compute statistical quantities and to 
compare the results to experiment. The results are shown to agree favorably with experi- 
ment. Quantitative agreement with particular experiments should not be expected due to 
differences in initial conditions such as shape and aspect ratio of the jet nozzle, intensity 
level and spectral content of the jet exit boundary layer and Reynolds number. The present 
results show significant influence of the spectral content of perturbation in the inlet mixing 
layer. Indeed, features of the jet flow such as the axis switching location and the length of 
the potential core vary by as much as 60% from experiment to experiment. 

Specific conclusions from the current study are outlined below. The first section 
describes conclusions about the numerical formulation, while the second section discusses 
conclusions based on the numerical simulation of complex jets. Contributions which make 
the present work unique are also highlighted. The final section outlines the suggestions for 
future work. 
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8.1 Numerical Formulation 


Numerical Approximation of Convection Terms 

It has been shown that many popular numerical approximations for convection terms 
lead to large dispersion and/or dissipation errors of the high-frequency modes, thus requir- 
ing upwards of 20-30 spatial points per wavelength for acceptable accuracy. Higher-order 
accurate compact schemes are used in this study which lead to greatly reduced dispersion 
errors and no dissipation errors. Solutions to relevant benchmark problems indicate that 6 
- 8 points per wavelength are sufficient with the fourth-order compact scheme for accept- 
able accuracy. The implicit treatment of derivatives results in a global formulation closer 
to spectral methods than explicit finite difference schemes. 

Enforcement of Continuity Equation with Higher-Order Compact Schemes 

While compact schemes have been widely used for compressible simulations, their use 
for incompressible flows is complicated by the lack of an evolution equation for pressure. 
A numerical formulation was developed in the present study which solves a Poisson equa- 
tion for pressure using a higher-order compact scheme. The improved accuracy and reso- 
lution characteristics was demonstrated through the solution of benchmark problems 
governed by the Euler and Navier-Stokes equations. It is confirmed that overall accuracy is 
limited by the weakest link. 


Extension to Curvilinear Grids 

The numerical formulation was extended to curvilinear grids which was used to con- 
centrate grid points in the mixing layers of the rectangular jet near the inflow. Downstream 
of the potential core, the flow is fully turbulent. Clustering of the grid to resolve large gra- 
dients due to small scale structures is not possible without the use of time-varying, adap- 
tive grids (which is outside the scope of the present study). Therefore, the grid is gradually 
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relaxed to uniform spacing in this region. Uniform grids enable manipulation of the dis- 
crete Poisson equation leading to a reduced stencil size and cost as compared with the cur- 
vilinear grids making the computational rate roughly one order of magnitude less. Part of 
the reason for the increase in computational rate is that convergence rates for the multigrid 
solution of the Poisson equation deteriorate with grid clustering and large aspect ratio. 
Therefore, the uniform grid formulation was utilized for most of the simulations presented 
in Chaps. 6 and 7. 

8.2 Simulation of Complex Jets 

Conclusions from the direct and large eddy simulations of complex jets are outlined in 
this section. Aspects of the current work which are different from the current state of the 
art are highlighted. 

Effect of Reynolds Number 

The results from the DNS of non-circular jets at low Reynolds numbers (Re = 750) in 
Chap. 6 show that many features of experiments at moderate to high Reynolds number are 
not captured. For example, low Reynolds number simulations in the present study and 
those by Miller et al (1995) show that the end of the potential core is not predicted within 
10 diameters, while experiments and the simulations at higher Reynolds number of the 
current study show an end at roughly x/D e - 4. Also transition to turbulent flow does not 
occur within the computational domain resulting in symmetrical structures and spikes in 
the velocity power spectra that are characteristic of laminar/transitional flow. It is clear 
that simulations at higher Reynolds numbers must be performed for comparison with 
experiment in the literature. 
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Streamwise Extent 


Current state of the art for the large eddy simulation of square and rectangular jets at 
moderate to high Reynolds number is that due to Grinstein and DeVore (1992), Grinstein 
(1993), and Grinstein (1996). The streamwise extent of those simulations covers the very 
near field (i.e., only the potential core region, x/D e < 5). The current study is unique in the 
sense that the near and medium fields are simulated. This allows the so-called characteris- 
tic decay region downstream of the potential core to be studied. In addition, the experi- 
mentally observed axis switching of the 2: 1 rectangular jet at x/D e = 7 is captured. 


Subgrid Scale Model 

The current study employs an explicit Smagorinsky subgrid scale model which allows 
the model effects to be quantified. Moreover, the effect of the unresolved scales was com- 
pared to those of the resolved scales by computing budget terms of the mean momentum 
equations. This analysis shows that the contribution of the unresolved stresses which are 
modeled are 15 times smaller than the those of the resolved stresses which are directly 
computed in the simulation. 

Budgets for Mean Momentum and Reynolds Stress Equations 

A numerical database of the time-dependent results from the simulations has been 
archived. The database was post-processed to compute statistical quantities and detailed 
budgets of the resolved mean momentum and Reynolds stress equations. In general, bal- 
ance of the terms of the mean momentum equation are quite good. Balance of terms in the 
Reynolds stress equations was incomplete because of the under-prediction of turbulent 
dissipation which occurs mostly at small scales not resolved in the simulation. The data- 
base can be used to aid in the turbulence modeling of complex jets. 
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Effect of Initial Conditions 

The results of the current study show that the initial conditions of the jet boundary 
layer at the inflow strongly effect the resulting dynamics. Forcing functions consisting of a 
single discrete mode show that the mixing layers in the potential core rolls up into planar 
vortex rings. At low Reynolds number, the fundamental mode saturates and decays, while 
at higher Reynolds number, transition to turbulence takes place. Axis-switching was 
observed in both cases. The addition of a subharmonic to the forcing function at low Rey- 
nolds number is shown to result in the formation of intense vortex ribs near the jet center- 
line which lead to a partial bifurcation of the jet. Axis switching was not observed in this 

case. 

Forcing functions based on random, broad band modes show that the roll up of the 
mixing layers in the potential core is fundamentally different than that due to discrete forc- 
ing. Planar vortex rings are not formed with broad mode forcing and a staggered roll up of 
the upper and lower mixing layer in the minor axis plane only is observed. As a result, 
mixing in the major axis plane is suppressed leading to switching of the major and minor 
axis at x/D e = 6.6 for the rectangular jet and x/D e = 5.9 for the elliptic jet at higher Rey- 
nolds number. Naturally developing non-circular jets (i.e. without strong discrete forcing) 
are better modeled with broad mode forcing. Therefore, simulations forced with a single 
sinusoidal mode leading to the formation of planar vortex rings provide an incomplete 
description of the axis switching phenomenon in naturally developing jets. The results of 
the present study provide such a picture. 


188 



Effect of Non-Uniform Boundary Layer Curvature 

Effects of non-uniform azimuthal curvature of the jet boundary layer at the nozzle are 
studied though large eddy simulations of rectangular, elliptic, and circular jets. The results 
show that the elliptic jet entrains more fluid into the jet core than the rectangular jet. 

Grid Resolution 

A grid resolution study was performed through simulations on 4 grid levels, by dou- 
bling the number of grid points on successive levels. The coarsest grid resolution (with 2.5 
points per streamwise fundamental wavelength) does not capture unsteadiness ( u ^ = 
0.005 at x/D e = 10) which is essential for accurate prediction of the potential core length 
and subsequent axis switching. Increasing the grid resolution to the second level results in 
an increased level of unsteadiness (u ms = 0.05 at x/D e = 10) but does not predict the end of 
the potential core or axis switching within the computational domain. The results from the 
level 3 grid (with 10 points per streamwise fundamental wavelength) show a well defined 
end of the potential core and axis switching in good agreement with experiment. The 
effect of increasing the grid resolution in the near field to the fourth level, does not signifi- 
cantly change those predictions. The main effect is that the prediction of the end of the 
potential core is shifted upstream by 0.7 diameters, which is most likely due to differences 
in the broad mode forcing function on the third and fourth grid levels. When the results on 
the fourth level are corrected for this difference, relatively grid independent results are 
achieved. The resolution requirement of roughly 10 streamwise grid points per wavelength 
reached through this exercise is consistent with conclusions from the solution of bench- 
mark problems using compact schemes in Chap. 5. 
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8.3 Future Work 


Future research directions are identified in this section. In the area of the numerical 
formulation, the efficiency of the Poisson solver with curvilinear grids could be improved 
which would lead to competitive computational rates as compared with the uniform grid 
formulation. This would be used take advantage of grid clustering capabilities in the very 
near field. The resolution after the end of the potential core would still be necessarily uni- 
form due to the presence of small scale turbulent structures. Time- varying, adaptive 
meshes would be required to track those structures. 

While the present numerical simulation extends the state of the art past the end of the 
potential core into the characteristic decay region, it would be desirable to simulate all 3 
regions of the non-circular jet; (i) potential core, (ii) characteristic decay, and (iii) axisym- 
metric region. It would also be desirable to increase resolution to the fourth grid level of 
the resolution study for the entire domain length. Code improvements would be required 
to reduce the required wall clock time and memory for the simulation such as the use of 
parallel processing and Fortran buffer in/ buffer out statements (reading and writing to 
disk to save memory allocation). Parallelization should be a major goal of future work. 
The simulations can also be extended to solve compressible flows for acoustic predictions. 
There, accurate solutions are required for the unsteady flow field which is then used as the 
source term of an acoustic analogy equation. 
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